The goal of the analyses in this notebook is to identify predictors of frog survival following translocation. The analyses are conducted in a Bayesian framework using the R package rstanarm.

Load packages

library(tidyverse)
Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ─────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.1 ──
✔ ggplot2 3.3.6     ✔ purrr   0.3.4
✔ tibble  3.1.7     ✔ dplyr   1.0.9
✔ tidyr   1.2.0     ✔ stringr 1.4.0
✔ readr   2.1.2     ✔ forcats 0.5.1
── Conflicts ────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
library(brms)
Loading required package: Rcpp
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio
Loading 'brms' package (version 2.17.0). Useful instructions
can be found by typing help('brms'). A more detailed introduction
to the package is available through vignette('brms_overview').

Attaching package: ‘brms’

The following object is masked from ‘package:stats’:

    ar
library(rstanarm)
This is rstanarm version 2.21.3
- See https://mc-stan.org/rstanarm/articles/priors for changes to default priors!
- Default priors may change, so it's safest to specify priors, even if equivalent to the defaults.
- For execution on a local, multicore CPU with excess RAM we recommend calling
  options(mc.cores = parallel::detectCores())

Attaching package: ‘rstanarm’

The following objects are masked from ‘package:brms’:

    dirichlet, exponential, get_y, lasso, ngrps
library(bayesplot)
This is bayesplot version 1.9.0
- Online documentation and vignettes at mc-stan.org/bayesplot
- bayesplot theme set to bayesplot::theme_default()
   * Does _not_ affect other ggplot2 plots
   * See ?bayesplot_theme_set for details on theme setting
library(tidybayes)

Attaching package: ‘tidybayes’

The following objects are masked from ‘package:brms’:

    dstudent_t, pstudent_t, qstudent_t, rstudent_t
library(loo)
This is loo version 2.5.1
- Online documentation and vignettes at mc-stan.org/loo
- As of v2.0.0 loo defaults to 1 core but we recommend using as many as possible. Use the 'cores' argument or set options(mc.cores = NUM_CORES) for an entire session. 
library(broom.mixed)
library(janitor)

Attaching package: ‘janitor’

The following objects are masked from ‘package:stats’:

    chisq.test, fisher.test
source("classification_summary.R")

Read in frog translocation dataset

d1 <- read_csv(here::here("data", "clean", "frog_translocation_final.csv")) %>% 
  drop_na(bd_load) %>%  # drop xx records from 70xxx where bd_load is NA
  mutate(
    first = if_else(order == 1, 1, 0),
    surv = if_else(survival < 0.5, 0, 1),
    male = if_else(sex == "m", 1, 0),
    elev_lo = if_else(elevation < 3020, 1, 0), # using interval defined by cut_interval, n = 2
    elev_z = arm::rescale(elevation),
    snowt_z = arm::rescale(snow_t),
    snowt1_z = arm::rescale(snow_t1),
    day_z = arm::rescale(day),
    length_z = arm::rescale(length),
    bdload_l = log10(bd_load + 1),
    bdload_z = arm::rescale(bdload_l),
    shore_c = shore - mean(shore),
    male_c = male - mean(male),
    elevlo_c = elev_lo - mean(elev_lo),
    first_c = first - mean(first),
    across(c(shore_c, male_c, elevlo_c, first_c), round, 3), 
    across(c(site_id, shore_c, first_c, donor, male_c, elevlo_c, year), as.factor),
    surv = as.integer(surv))
Rows: 779 Columns: 16── Column specification ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr   (2): pit_tag_ref, sex
dbl  (13): site_id, elevation, shore, snow_t, snow_t1, year, day, order, donor, length, weight, bd_load, survival
date  (1): date
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# Add a descriptive translocation_id
d1 <- d1 %>% 
  distinct(site_id, year) %>% 
  arrange(site_id, year) %>% 
  group_by(site_id) %>% 
  mutate(transno = seq_len(n())) %>% 
  ungroup(site_id) %>% 
  unite(translocation_id, c("site_id", "transno"), remove = FALSE) %>% 
  select(-transno) %>% 
  mutate(translocation_id = as.factor(translocation_id)) %>% 
  inner_join(d1, by = c("site_id", "year")) %>% 
  relocate(translocation_id, .after = site_id)

Null values - check

d1 %>% summarize(across(everything(), ~ sum(is.na(.))))

Model using non-standardized predictors

Model specification

m1a <- stan_glm(
  surv ~ length + snow_t + snow_t1 + day + bdload_l + elevation + first + male + donor,
  data = d1,
  family = binomial,
  chains = 4, 
  iter = 5000*2,
  cores = 4,
  seed = 84735)
starting worker pid=10413 on localhost:11976 at 13:34:41.032
starting worker pid=10442 on localhost:11976 at 13:34:41.153
starting worker pid=10471 on localhost:11976 at 13:34:41.276
starting worker pid=10500 on localhost:11976 at 13:34:41.395

SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 7.5e-05 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.75 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 10000 [  0%]  (Warmup)
Chain 1: Iteration: 1000 / 10000 [ 10%]  (Warmup)

SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 2.9e-05 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.29 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 10000 [  0%]  (Warmup)
Chain 1: Iteration: 2000 / 10000 [ 20%]  (Warmup)
Chain 2: Iteration: 1000 / 10000 [ 10%]  (Warmup)

SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 3.1e-05 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.31 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:    1 / 10000 [  0%]  (Warmup)
Chain 1: Iteration: 3000 / 10000 [ 30%]  (Warmup)
Chain 2: Iteration: 2000 / 10000 [ 20%]  (Warmup)

SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 2.8e-05 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.28 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:    1 / 10000 [  0%]  (Warmup)
Chain 3: Iteration: 1000 / 10000 [ 10%]  (Warmup)
Chain 1: Iteration: 4000 / 10000 [ 40%]  (Warmup)
Chain 2: Iteration: 3000 / 10000 [ 30%]  (Warmup)
Chain 4: Iteration: 1000 / 10000 [ 10%]  (Warmup)
Chain 3: Iteration: 2000 / 10000 [ 20%]  (Warmup)
Chain 1: Iteration: 5000 / 10000 [ 50%]  (Warmup)
Chain 1: Iteration: 5001 / 10000 [ 50%]  (Sampling)
Chain 2: Iteration: 4000 / 10000 [ 40%]  (Warmup)
Chain 4: Iteration: 2000 / 10000 [ 20%]  (Warmup)
Chain 3: Iteration: 3000 / 10000 [ 30%]  (Warmup)
Chain 1: Iteration: 6000 / 10000 [ 60%]  (Sampling)
Chain 2: Iteration: 5000 / 10000 [ 50%]  (Warmup)
Chain 2: Iteration: 5001 / 10000 [ 50%]  (Sampling)
Chain 4: Iteration: 3000 / 10000 [ 30%]  (Warmup)
Chain 3: Iteration: 4000 / 10000 [ 40%]  (Warmup)
Chain 1: Iteration: 7000 / 10000 [ 70%]  (Sampling)
Chain 2: Iteration: 6000 / 10000 [ 60%]  (Sampling)
Chain 4: Iteration: 4000 / 10000 [ 40%]  (Warmup)
Chain 3: Iteration: 5000 / 10000 [ 50%]  (Warmup)
Chain 3: Iteration: 5001 / 10000 [ 50%]  (Sampling)
Chain 1: Iteration: 8000 / 10000 [ 80%]  (Sampling)
Chain 2: Iteration: 7000 / 10000 [ 70%]  (Sampling)
Chain 4: Iteration: 5000 / 10000 [ 50%]  (Warmup)
Chain 4: Iteration: 5001 / 10000 [ 50%]  (Sampling)
Chain 3: Iteration: 6000 / 10000 [ 60%]  (Sampling)
Chain 1: Iteration: 9000 / 10000 [ 90%]  (Sampling)
Chain 2: Iteration: 8000 / 10000 [ 80%]  (Sampling)
Chain 4: Iteration: 6000 / 10000 [ 60%]  (Sampling)
Chain 3: Iteration: 7000 / 10000 [ 70%]  (Sampling)
Chain 1: Iteration: 10000 / 10000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 1.04962 seconds (Warm-up)
Chain 1:                1.09124 seconds (Sampling)
Chain 1:                2.14086 seconds (Total)
Chain 1: 
Chain 2: Iteration: 9000 / 10000 [ 90%]  (Sampling)
Chain 4: Iteration: 7000 / 10000 [ 70%]  (Sampling)
Chain 3: Iteration: 8000 / 10000 [ 80%]  (Sampling)
Chain 2: Iteration: 10000 / 10000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 1.01617 seconds (Warm-up)
Chain 2:                1.10254 seconds (Sampling)
Chain 2:                2.11871 seconds (Total)
Chain 2: 
Chain 4: Iteration: 8000 / 10000 [ 80%]  (Sampling)
Chain 3: Iteration: 9000 / 10000 [ 90%]  (Sampling)
Chain 4: Iteration: 9000 / 10000 [ 90%]  (Sampling)
Chain 3: Iteration: 10000 / 10000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 1.01876 seconds (Warm-up)
Chain 3:                1.17841 seconds (Sampling)
Chain 3:                2.19717 seconds (Total)
Chain 3: 
Chain 4: Iteration: 10000 / 10000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 1.03517 seconds (Warm-up)
Chain 4:                1.15732 seconds (Sampling)
Chain 4:                2.1925 seconds (Total)
Chain 4: 

Priors & chain diagnostics

prior_summary(m1a)
Priors for model 'm1a' 
------
Intercept (after predictors centered)
 ~ normal(location = 0, scale = 2.5)

Coefficients
  Specified prior:
    ~ normal(location = [0,0,0,...], scale = [2.5,2.5,2.5,...])
  Adjusted prior:
    ~ normal(location = [0,0,0,...], scale = [0.526,0.049,0.073,...])
------
See help('prior_summary.stanreg') for more details
mcmc_trace(m1a, size = 0.1)

mcmc_dens_overlay(m1a)


summary(m1a)

Model Info:
 function:     stan_glm
 family:       binomial [logit]
 formula:      surv ~ length + snow_t + snow_t1 + day + bdload_l + elevation + 
       first + male + donor
 algorithm:    sampling
 sample:       20000 (posterior sample size)
 priors:       see help('prior_summary')
 observations: 767
 predictors:   11

Estimates:
              mean   sd    10%   50%   90%
(Intercept)  -8.5    2.3 -11.5  -8.5  -5.5
length        0.0    0.0   0.0   0.0   0.0
snow_t        0.0    0.0   0.0   0.0   0.0
snow_t1       0.0    0.0   0.0   0.0   0.0
day           0.0    0.0   0.0   0.0   0.0
bdload_l     -0.1    0.1  -0.1  -0.1   0.0
elevation     0.0    0.0   0.0   0.0   0.0
first         0.4    0.2   0.1   0.4   0.7
male          0.3    0.2   0.1   0.3   0.5
donor70567   -1.3    0.4  -1.9  -1.3  -0.8
donor72996   -2.2    0.4  -2.7  -2.2  -1.7

Fit Diagnostics:
           mean   sd   10%   50%   90%
mean_PPD 0.3    0.0  0.3   0.3   0.3  

The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).

MCMC diagnostics
              mcse Rhat n_eff
(Intercept)   0.0  1.0  24916
length        0.0  1.0  25653
snow_t        0.0  1.0  22266
snow_t1       0.0  1.0  21140
day           0.0  1.0  22645
bdload_l      0.0  1.0  29642
elevation     0.0  1.0  24654
first         0.0  1.0  23563
male          0.0  1.0  28297
donor70567    0.0  1.0  20512
donor72996    0.0  1.0  18851
mean_PPD      0.0  1.0  23850
log-posterior 0.0  1.0   8254

For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).
  • mcmc diagnostics all look good.

List model variable names

get_variables(m1a)
 [1] "(Intercept)"   "length"        "snow_t"        "snow_t1"       "day"           "bdload_l"      "elevation"     "first"        
 [9] "male"          "donor70567"    "donor72996"    "accept_stat__" "stepsize__"    "treedepth__"   "n_leapfrog__"  "divergent__"  
[17] "energy__"     

Posterior predictive check

pp_check(m1a, plotfun = "stat") +
    xlab("Frog survival rate")

  • Posterior predictive distribution implied by regression model closely matches original data (based on mean frog survival).

Calculate leave-one-out cross-validation (loo)

loo_m1a <- loo(m1a, cores = 4, save_psis = TRUE)
print(loo_m1a)

Computed from 20000 by 767 log-likelihood matrix
------
Monte Carlo SE of elpd_loo is 0.0.

All Pareto k estimates are good (k < 0.5).
See help('pareto-k-diagnostic') for details.

Posterior analysis

tidy(m1a, conf.int = TRUE, conf.level = 0.90)

mcmc_intervals(m1a, point_est = "median", prob = 0.80, prob_outer = 0.9) +
  geom_vline(xintercept = 0, color = "red", linetype = "dashed", lwd = 0.5) +  # add 0-line to accentuate default line
  ggtitle("Model with non-standardized predictors only",
          "Posterior distributions with medians & 90% uncertainty intervals")


mcmc_areas(m1a, area_method = "equal height", prob = 0.90) +
  geom_vline(xintercept = 0, color = "red", linetype = "dashed", lwd = 0.5) +  
  ggtitle("Model with non-standardized predictors only",
          "Posterior distributions with medians & 90% uncertainty intervals")


mean(bayes_R2(m1a))
[1] 0.2728229
  • Using non-standardized predictors makes plot hard to read due to vastly different scales of predictors.
  • Based on 90% uncertainty intervals, important predictors are snow_t1, elevation, first, donor70567, and donor72996.

Model using standardized predictors

Model specification

m1b <- stan_glm(
  surv ~ length_z + snowt_z + snowt1_z + day_z + bdload_z + elev_z + first_c + male_c + donor,
  data = d1,
  family = binomial,
  chains = 4, 
  iter = 5000*2, 
  cores = 4, 
  seed = 84735)
starting worker pid=10862 on localhost:11976 at 13:35:10.488
starting worker pid=10891 on localhost:11976 at 13:35:10.614
starting worker pid=10920 on localhost:11976 at 13:35:10.733
starting worker pid=10949 on localhost:11976 at 13:35:10.863

SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 4.8e-05 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.48 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 10000 [  0%]  (Warmup)
Chain 1: Iteration: 1000 / 10000 [ 10%]  (Warmup)

SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 2.7e-05 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.27 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 10000 [  0%]  (Warmup)
Chain 1: Iteration: 2000 / 10000 [ 20%]  (Warmup)
Chain 2: Iteration: 1000 / 10000 [ 10%]  (Warmup)

SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 2.6e-05 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.26 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:    1 / 10000 [  0%]  (Warmup)
Chain 1: Iteration: 3000 / 10000 [ 30%]  (Warmup)
Chain 2: Iteration: 2000 / 10000 [ 20%]  (Warmup)

SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 3.2e-05 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.32 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:    1 / 10000 [  0%]  (Warmup)
Chain 3: Iteration: 1000 / 10000 [ 10%]  (Warmup)
Chain 1: Iteration: 4000 / 10000 [ 40%]  (Warmup)
Chain 2: Iteration: 3000 / 10000 [ 30%]  (Warmup)
Chain 4: Iteration: 1000 / 10000 [ 10%]  (Warmup)
Chain 3: Iteration: 2000 / 10000 [ 20%]  (Warmup)
Chain 1: Iteration: 5000 / 10000 [ 50%]  (Warmup)
Chain 1: Iteration: 5001 / 10000 [ 50%]  (Sampling)
Chain 2: Iteration: 4000 / 10000 [ 40%]  (Warmup)
Chain 4: Iteration: 2000 / 10000 [ 20%]  (Warmup)
Chain 3: Iteration: 3000 / 10000 [ 30%]  (Warmup)
Chain 2: Iteration: 5000 / 10000 [ 50%]  (Warmup)
Chain 2: Iteration: 5001 / 10000 [ 50%]  (Sampling)
Chain 1: Iteration: 6000 / 10000 [ 60%]  (Sampling)
Chain 4: Iteration: 3000 / 10000 [ 30%]  (Warmup)
Chain 3: Iteration: 4000 / 10000 [ 40%]  (Warmup)
Chain 2: Iteration: 6000 / 10000 [ 60%]  (Sampling)
Chain 1: Iteration: 7000 / 10000 [ 70%]  (Sampling)
Chain 4: Iteration: 4000 / 10000 [ 40%]  (Warmup)
Chain 3: Iteration: 5000 / 10000 [ 50%]  (Warmup)
Chain 3: Iteration: 5001 / 10000 [ 50%]  (Sampling)
Chain 2: Iteration: 7000 / 10000 [ 70%]  (Sampling)
Chain 4: Iteration: 5000 / 10000 [ 50%]  (Warmup)
Chain 4: Iteration: 5001 / 10000 [ 50%]  (Sampling)
Chain 1: Iteration: 8000 / 10000 [ 80%]  (Sampling)
Chain 3: Iteration: 6000 / 10000 [ 60%]  (Sampling)
Chain 2: Iteration: 8000 / 10000 [ 80%]  (Sampling)
Chain 1: Iteration: 9000 / 10000 [ 90%]  (Sampling)
Chain 4: Iteration: 6000 / 10000 [ 60%]  (Sampling)
Chain 3: Iteration: 7000 / 10000 [ 70%]  (Sampling)
Chain 2: Iteration: 9000 / 10000 [ 90%]  (Sampling)
Chain 1: Iteration: 10000 / 10000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 1.04687 seconds (Warm-up)
Chain 1:                1.09385 seconds (Sampling)
Chain 1:                2.14072 seconds (Total)
Chain 1: 
Chain 4: Iteration: 7000 / 10000 [ 70%]  (Sampling)
Chain 3: Iteration: 8000 / 10000 [ 80%]  (Sampling)
Chain 2: Iteration: 10000 / 10000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 1.0252 seconds (Warm-up)
Chain 2:                1.10545 seconds (Sampling)
Chain 2:                2.13065 seconds (Total)
Chain 2: 
Chain 4: Iteration: 8000 / 10000 [ 80%]  (Sampling)
Chain 3: Iteration: 9000 / 10000 [ 90%]  (Sampling)
Chain 4: Iteration: 9000 / 10000 [ 90%]  (Sampling)
Chain 3: Iteration: 10000 / 10000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 1.04144 seconds (Warm-up)
Chain 3:                1.18587 seconds (Sampling)
Chain 3:                2.22731 seconds (Total)
Chain 3: 
Chain 4: Iteration: 10000 / 10000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 1.04213 seconds (Warm-up)
Chain 4:                1.14531 seconds (Sampling)
Chain 4:                2.18744 seconds (Total)
Chain 4: 

Priors, chain diagnostics, posterior predictive check

prior_summary(m1b)
Priors for model 'm1b' 
------
Intercept (after predictors centered)
 ~ normal(location = 0, scale = 2.5)

Coefficients
  Specified prior:
    ~ normal(location = [0,0,0,...], scale = [2.5,2.5,2.5,...])
  Adjusted prior:
    ~ normal(location = [0,0,0,...], scale = [5.00,5.00,5.00,...])
------
See help('prior_summary.stanreg') for more details
mcmc_trace(m1b, size = 0.1)

mcmc_dens_overlay(m1a)


summary(m1b)

Model Info:
 function:     stan_glm
 family:       binomial [logit]
 formula:      surv ~ length_z + snowt_z + snowt1_z + day_z + bdload_z + elev_z + 
       first_c + male_c + donor
 algorithm:    sampling
 sample:       20000 (posterior sample size)
 priors:       see help('prior_summary')
 observations: 767
 predictors:   11

Estimates:
               mean   sd   10%   50%   90%
(Intercept)   0.4    0.4  0.0   0.4   0.9 
length_z      0.0    0.2 -0.3   0.0   0.3 
snowt_z       0.0    0.2 -0.3   0.0   0.3 
snowt1_z     -1.3    0.2 -1.6  -1.3  -1.0 
day_z         0.2    0.2 -0.1   0.2   0.4 
bdload_z     -0.2    0.2 -0.4  -0.2   0.1 
elev_z        1.5    0.3  1.2   1.5   1.9 
first_c0.353  0.4    0.2  0.1   0.4   0.7 
male_c0.548   0.3    0.2  0.1   0.3   0.5 
donor70567   -1.3    0.4 -1.9  -1.3  -0.8 
donor72996   -2.2    0.4 -2.7  -2.2  -1.7 

Fit Diagnostics:
           mean   sd   10%   50%   90%
mean_PPD 0.3    0.0  0.3   0.3   0.3  

The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).

MCMC diagnostics
              mcse Rhat n_eff
(Intercept)   0.0  1.0  21959
length_z      0.0  1.0  25653
snowt_z       0.0  1.0  22266
snowt1_z      0.0  1.0  21140
day_z         0.0  1.0  22645
bdload_z      0.0  1.0  29642
elev_z        0.0  1.0  24654
first_c0.353  0.0  1.0  23563
male_c0.548   0.0  1.0  28297
donor70567    0.0  1.0  20512
donor72996    0.0  1.0  18851
mean_PPD      0.0  1.0  23850
log-posterior 0.0  1.0   8254

For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).
  • mcmc diagnostics all look good.

List model variable names

get_variables(m1b)
 [1] "(Intercept)"   "length_z"      "snowt_z"       "snowt1_z"      "day_z"         "bdload_z"      "elev_z"        "first_c0.353" 
 [9] "male_c0.548"   "donor70567"    "donor72996"    "accept_stat__" "stepsize__"    "treedepth__"   "n_leapfrog__"  "divergent__"  
[17] "energy__"     

Posterior predictive check

pp_check(m1b, plotfun = "stat") +
    xlab("Frog survival rate")

  • Posterior predictive distribution implied by regression model closely matches original data (based on mean frog survival).

Calculate loo

loo_m1b <- loo(m1b, cores = 4, save_psis = TRUE)
print(loo_m1b)

Computed from 20000 by 767 log-likelihood matrix
------
Monte Carlo SE of elpd_loo is 0.0.

All Pareto k estimates are good (k < 0.5).
See help('pareto-k-diagnostic') for details.

Posterior analysis

tidy(m1b, conf.int = TRUE, conf.level = 0.90)

mcmc_areas(m1b, area_method = "equal height", prob = 0.90) +
  geom_vline(xintercept = 0, color = "red", linetype = "dashed", lwd = 0.5) +  
  ggtitle("Model with standardized predictors only",
          "Posterior distributions with medians & 90% uncertainty intervals")


mean(bayes_R2(m1b))
[1] 0.2728229
  • As expected, using standardized predictors allows easy comparison across predictors. Important predictors are unchanged.
  • Based on 90% uncertainty intervals, important predictors are snow_t1, elevation, first, donor70567, and donor72996.

Compare models

loo_compare(loo_m1a, loo_m1b)
    elpd_diff se_diff
m1a 0.0       0.0    
m1b 0.0       0.0    

Possible grouping structures to add to model

The dataset has three groups: site_id, translocation_id, and donor site. Donor site is not suitable as a grouping variable because there are only 3 levels (i.e., site ids), not enough to make this useful. Site_id and translocation_id contain 12 and 24 levels, respectively, making them suitable as grouping variables. Including both site_id and translocation_id as nested grouping variables (i.e., translocations nested within site ids) also seems justified. Having more levels within a group would seem to increase the power to detect effects of predictors. With fewer levels/less variability, effect may be accounted for primarily by grouping variable itself instead of predictor. With only 1-4 translocations per site and only 12 site ids total, it is possible that model with site_id as sole grouping variable and model with site_id and translocation_id as nested grouping variables will attribute more of the variation in frog survival to the grouping variables and less to the predictors. A priori, that would suggest that the model with translocation_id as the grouping variable might be the most useful model structure (assuming that all models have similar posterior predictive accuracy). Lacking any strong rationale for one model structure over another, built 3 models and compared their posterior predictive accuracy using loo.

Model using standardized predictors and site_id grouping variable

Model specification

m1c <- stan_glmer(
  surv ~ length_z + snowt_z + snowt1_z + day_z + bdload_z + elev_z + first_c + male_c + donor + (1 | site_id),
  data = d1,
  family = binomial,
  chains = 4, 
  iter = 5000*2, 
  cores = 4,
  seed = 84735)
starting worker pid=11308 on localhost:11976 at 13:35:38.500
starting worker pid=11337 on localhost:11976 at 13:35:38.651
starting worker pid=11372 on localhost:11976 at 13:35:38.790
starting worker pid=11403 on localhost:11976 at 13:35:38.937

SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 1).
Chain 1: Rejecting initial value:
Chain 1:   Log probability evaluates to log(0), i.e. negative infinity.
Chain 1:   Stan can't start sampling from this initial value.
Chain 1: 
Chain 1: Gradient evaluation took 0.000169 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.69 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 10000 [  0%]  (Warmup)

SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 0.000109 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 1.09 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 10000 [  0%]  (Warmup)

SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 0.000118 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 1.18 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:    1 / 10000 [  0%]  (Warmup)

SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 4).
Chain 4: Rejecting initial value:
Chain 4:   Log probability evaluates to log(0), i.e. negative infinity.
Chain 4:   Stan can't start sampling from this initial value.
Chain 4: 
Chain 4: Gradient evaluation took 0.000119 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 1.19 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:    1 / 10000 [  0%]  (Warmup)
Chain 3: Iteration: 1000 / 10000 [ 10%]  (Warmup)
Chain 2: Iteration: 1000 / 10000 [ 10%]  (Warmup)
Chain 1: Iteration: 1000 / 10000 [ 10%]  (Warmup)
Chain 4: Iteration: 1000 / 10000 [ 10%]  (Warmup)
Chain 3: Iteration: 2000 / 10000 [ 20%]  (Warmup)
Chain 1: Iteration: 2000 / 10000 [ 20%]  (Warmup)
Chain 2: Iteration: 2000 / 10000 [ 20%]  (Warmup)
Chain 4: Iteration: 2000 / 10000 [ 20%]  (Warmup)
Chain 3: Iteration: 3000 / 10000 [ 30%]  (Warmup)
Chain 1: Iteration: 3000 / 10000 [ 30%]  (Warmup)
Chain 2: Iteration: 3000 / 10000 [ 30%]  (Warmup)
Chain 4: Iteration: 3000 / 10000 [ 30%]  (Warmup)
Chain 3: Iteration: 4000 / 10000 [ 40%]  (Warmup)
Chain 1: Iteration: 4000 / 10000 [ 40%]  (Warmup)
Chain 2: Iteration: 4000 / 10000 [ 40%]  (Warmup)
Chain 4: Iteration: 4000 / 10000 [ 40%]  (Warmup)
Chain 2: Iteration: 5000 / 10000 [ 50%]  (Warmup)
Chain 2: Iteration: 5001 / 10000 [ 50%]  (Sampling)
Chain 1: Iteration: 5000 / 10000 [ 50%]  (Warmup)
Chain 1: Iteration: 5001 / 10000 [ 50%]  (Sampling)
Chain 3: Iteration: 5000 / 10000 [ 50%]  (Warmup)
Chain 3: Iteration: 5001 / 10000 [ 50%]  (Sampling)
Chain 4: Iteration: 5000 / 10000 [ 50%]  (Warmup)
Chain 4: Iteration: 5001 / 10000 [ 50%]  (Sampling)
Chain 2: Iteration: 6000 / 10000 [ 60%]  (Sampling)
Chain 3: Iteration: 6000 / 10000 [ 60%]  (Sampling)
Chain 1: Iteration: 6000 / 10000 [ 60%]  (Sampling)
Chain 4: Iteration: 6000 / 10000 [ 60%]  (Sampling)
Chain 2: Iteration: 7000 / 10000 [ 70%]  (Sampling)
Chain 3: Iteration: 7000 / 10000 [ 70%]  (Sampling)
Chain 4: Iteration: 7000 / 10000 [ 70%]  (Sampling)
Chain 1: Iteration: 7000 / 10000 [ 70%]  (Sampling)
Chain 2: Iteration: 8000 / 10000 [ 80%]  (Sampling)
Chain 3: Iteration: 8000 / 10000 [ 80%]  (Sampling)
Chain 4: Iteration: 8000 / 10000 [ 80%]  (Sampling)
Chain 1: Iteration: 8000 / 10000 [ 80%]  (Sampling)
Chain 2: Iteration: 9000 / 10000 [ 90%]  (Sampling)
Chain 4: Iteration: 9000 / 10000 [ 90%]  (Sampling)
Chain 3: Iteration: 9000 / 10000 [ 90%]  (Sampling)
Chain 1: Iteration: 9000 / 10000 [ 90%]  (Sampling)
Chain 2: Iteration: 10000 / 10000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 26.8153 seconds (Warm-up)
Chain 2:                25.8212 seconds (Sampling)
Chain 2:                52.6365 seconds (Total)
Chain 2: 
Chain 4: Iteration: 10000 / 10000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 27.3018 seconds (Warm-up)
Chain 4:                25.8839 seconds (Sampling)
Chain 4:                53.1858 seconds (Total)
Chain 4: 
Chain 3: Iteration: 10000 / 10000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 26.6674 seconds (Warm-up)
Chain 3:                27.4283 seconds (Sampling)
Chain 3:                54.0958 seconds (Total)
Chain 3: 
Chain 1: Iteration: 10000 / 10000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 27.1176 seconds (Warm-up)
Chain 1:                30.398 seconds (Sampling)
Chain 1:                57.5155 seconds (Total)
Chain 1: 

Priors & chain diagnostics

prior_summary(m1c)
Priors for model 'm1c' 
------
Intercept (after predictors centered)
 ~ normal(location = 0, scale = 2.5)

Coefficients
  Specified prior:
    ~ normal(location = [0,0,0,...], scale = [2.5,2.5,2.5,...])
  Adjusted prior:
    ~ normal(location = [0,0,0,...], scale = [5.00,5.00,5.00,...])

Covariance
 ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
------
See help('prior_summary.stanreg') for more details
b1 <- c("(Intercept)", "length_z", "snowt_z", "snowt1_z", "day_z", "bdload_z", "elev_z", "first_c0.353", "male_c0.548", "donor70567", "donor72996")

mcmc_trace(m1c, size = 0.1, pars = b1)

mcmc_dens_overlay(m1c, pars = b1)


summary(m1c)

Model Info:
 function:     stan_glmer
 family:       binomial [logit]
 formula:      surv ~ length_z + snowt_z + snowt1_z + day_z + bdload_z + elev_z + 
       first_c + male_c + donor + (1 | site_id)
 algorithm:    sampling
 sample:       20000 (posterior sample size)
 priors:       see help('prior_summary')
 observations: 767
 groups:       site_id (12)

Estimates:
                                         mean   sd   10%   50%   90%
(Intercept)                             0.0    1.2 -1.5   0.0   1.3 
length_z                                0.0    0.2 -0.3   0.0   0.3 
snowt_z                                 0.0    0.4 -0.5   0.0   0.5 
snowt1_z                                0.0    0.5 -0.7  -0.1   0.7 
day_z                                   0.2    0.4 -0.4   0.2   0.7 
bdload_z                               -0.2    0.2 -0.5  -0.2   0.1 
elev_z                                  1.7    1.0  0.5   1.7   3.0 
first_c0.353                            0.3    0.2  0.0   0.3   0.7 
male_c0.548                             0.4    0.2  0.1   0.3   0.6 
donor70567                             -0.8    1.7 -2.8  -0.8   1.3 
donor72996                             -1.9    1.3 -3.4  -1.9  -0.3 
b[(Intercept) site_id:70134]           -0.4    0.6 -1.2  -0.4   0.3 
b[(Intercept) site_id:70370]           -1.7    1.0 -2.9  -1.6  -0.5 
b[(Intercept) site_id:70413]            0.0    1.3 -1.6   0.0   1.5 
b[(Intercept) site_id:70414]           -1.2    1.2 -2.7  -1.0   0.2 
b[(Intercept) site_id:70449]            1.5    0.6  0.8   1.5   2.3 
b[(Intercept) site_id:70505]            0.2    0.6 -0.5   0.2   1.0 
b[(Intercept) site_id:70550]            0.6    0.6 -0.2   0.5   1.4 
b[(Intercept) site_id:70556]           -0.9    1.0 -2.2  -0.9   0.4 
b[(Intercept) site_id:70619]           -0.6    0.7 -1.5  -0.5   0.3 
b[(Intercept) site_id:70628]            1.1    0.8  0.2   1.1   2.1 
b[(Intercept) site_id:70641]            0.0    0.8 -0.9   0.0   1.0 
b[(Intercept) site_id:74976]            1.0    1.0 -0.2   0.9   2.3 
Sigma[site_id:(Intercept),(Intercept)]  1.9    1.4  0.6   1.5   3.5 

Fit Diagnostics:
           mean   sd   10%   50%   90%
mean_PPD 0.3    0.0  0.3   0.3   0.3  

The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).

MCMC diagnostics
                                       mcse Rhat n_eff
(Intercept)                            0.0  1.0   9604
length_z                               0.0  1.0  29640
snowt_z                                0.0  1.0  12471
snowt1_z                               0.0  1.0   9663
day_z                                  0.0  1.0   9895
bdload_z                               0.0  1.0  29229
elev_z                                 0.0  1.0   8370
first_c0.353                           0.0  1.0  25887
male_c0.548                            0.0  1.0  30649
donor70567                             0.0  1.0  12665
donor72996                             0.0  1.0   9205
b[(Intercept) site_id:70134]           0.0  1.0  12243
b[(Intercept) site_id:70370]           0.0  1.0   9052
b[(Intercept) site_id:70413]           0.0  1.0  16783
b[(Intercept) site_id:70414]           0.0  1.0  13633
b[(Intercept) site_id:70449]           0.0  1.0   9662
b[(Intercept) site_id:70505]           0.0  1.0  10089
b[(Intercept) site_id:70550]           0.0  1.0   9125
b[(Intercept) site_id:70556]           0.0  1.0  13480
b[(Intercept) site_id:70619]           0.0  1.0   9114
b[(Intercept) site_id:70628]           0.0  1.0   9365
b[(Intercept) site_id:70641]           0.0  1.0   8782
b[(Intercept) site_id:74976]           0.0  1.0  10720
Sigma[site_id:(Intercept),(Intercept)] 0.0  1.0   7013
mean_PPD                               0.0  1.0  18973
log-posterior                          0.1  1.0   5668

For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).

Posterior predictive check

pp_check(m1c, plotfun = "stat") +
    xlab("Frog survival rate")

Calculate loo

loo_m1c <- loo(m1c, cores = 4, save_psis = TRUE)
print(loo_m1c)

Computed from 20000 by 767 log-likelihood matrix
------
Monte Carlo SE of elpd_loo is 0.0.

All Pareto k estimates are good (k < 0.5).
See help('pareto-k-diagnostic') for details.

Posterior analysis

tidy(m1c, conf.int = TRUE, conf.level = 0.90)

mcmc_areas(m1c, area_method = "equal height", prob = 0.90, pars = b1) +
  geom_vline(xintercept = 0, color = "red", linetype = "dashed", lwd = 0.5) +  
  ggtitle("Model with standardized predictors & site_id as grouping variable",
          "Posterior distributions with medians & 90% uncertainty intervals")


mean(bayes_R2(m1c))
[1] 0.3169546
# Frog survival by site_id
mcmc_areas_ridges(m1c, regex_pars = "site_id") +
  geom_vline(xintercept = 0, color = "red", linetype = "dashed", lwd = 0.5) +  
  ggtitle("Frog survival for each recipient site_id")

  • Based on 90% uncertainty intervals, important predictors reduced to elevation and male.

Model using standardized predictors and translocation_id grouping variable)

Model specification

m1d <- stan_glmer(
  surv ~ length_z + snowt_z + snowt1_z + day_z + bdload_z + elev_z + first_c + male_c + donor + (1 | translocation_id),
  data = d1,
  family = binomial,
  chains = 4, 
  iter = 5000*2, 
  cores = 4,
  seed = 84735)
starting worker pid=11797 on localhost:11976 at 13:37:04.696
starting worker pid=11826 on localhost:11976 at 13:37:04.820
starting worker pid=11855 on localhost:11976 at 13:37:04.938
starting worker pid=11884 on localhost:11976 at 13:37:05.060

SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 1).
Chain 1: Rejecting initial value:
Chain 1:   Log probability evaluates to log(0), i.e. negative infinity.
Chain 1:   Stan can't start sampling from this initial value.
Chain 1: 
Chain 1: Gradient evaluation took 0.000111 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.11 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 10000 [  0%]  (Warmup)

SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 0.000108 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 1.08 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 10000 [  0%]  (Warmup)

SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 0.00012 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 1.2 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:    1 / 10000 [  0%]  (Warmup)

SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 0.000115 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 1.15 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:    1 / 10000 [  0%]  (Warmup)
Chain 2: Iteration: 1000 / 10000 [ 10%]  (Warmup)
Chain 1: Iteration: 1000 / 10000 [ 10%]  (Warmup)
Chain 3: Iteration: 1000 / 10000 [ 10%]  (Warmup)
Chain 4: Iteration: 1000 / 10000 [ 10%]  (Warmup)
Chain 2: Iteration: 2000 / 10000 [ 20%]  (Warmup)
Chain 1: Iteration: 2000 / 10000 [ 20%]  (Warmup)
Chain 3: Iteration: 2000 / 10000 [ 20%]  (Warmup)
Chain 4: Iteration: 2000 / 10000 [ 20%]  (Warmup)
Chain 2: Iteration: 3000 / 10000 [ 30%]  (Warmup)
Chain 1: Iteration: 3000 / 10000 [ 30%]  (Warmup)
Chain 3: Iteration: 3000 / 10000 [ 30%]  (Warmup)
Chain 4: Iteration: 3000 / 10000 [ 30%]  (Warmup)
Chain 2: Iteration: 4000 / 10000 [ 40%]  (Warmup)
Chain 1: Iteration: 4000 / 10000 [ 40%]  (Warmup)
Chain 3: Iteration: 4000 / 10000 [ 40%]  (Warmup)
Chain 4: Iteration: 4000 / 10000 [ 40%]  (Warmup)
Chain 2: Iteration: 5000 / 10000 [ 50%]  (Warmup)
Chain 2: Iteration: 5001 / 10000 [ 50%]  (Sampling)
Chain 1: Iteration: 5000 / 10000 [ 50%]  (Warmup)
Chain 1: Iteration: 5001 / 10000 [ 50%]  (Sampling)
Chain 3: Iteration: 5000 / 10000 [ 50%]  (Warmup)
Chain 3: Iteration: 5001 / 10000 [ 50%]  (Sampling)
Chain 4: Iteration: 5000 / 10000 [ 50%]  (Warmup)
Chain 4: Iteration: 5001 / 10000 [ 50%]  (Sampling)
Chain 2: Iteration: 6000 / 10000 [ 60%]  (Sampling)
Chain 1: Iteration: 6000 / 10000 [ 60%]  (Sampling)
Chain 3: Iteration: 6000 / 10000 [ 60%]  (Sampling)
Chain 4: Iteration: 6000 / 10000 [ 60%]  (Sampling)
Chain 2: Iteration: 7000 / 10000 [ 70%]  (Sampling)
Chain 1: Iteration: 7000 / 10000 [ 70%]  (Sampling)
Chain 3: Iteration: 7000 / 10000 [ 70%]  (Sampling)
Chain 4: Iteration: 7000 / 10000 [ 70%]  (Sampling)
Chain 2: Iteration: 8000 / 10000 [ 80%]  (Sampling)
Chain 1: Iteration: 8000 / 10000 [ 80%]  (Sampling)
Chain 3: Iteration: 8000 / 10000 [ 80%]  (Sampling)
Chain 4: Iteration: 8000 / 10000 [ 80%]  (Sampling)
Chain 1: Iteration: 9000 / 10000 [ 90%]  (Sampling)
Chain 2: Iteration: 9000 / 10000 [ 90%]  (Sampling)
Chain 3: Iteration: 9000 / 10000 [ 90%]  (Sampling)
Chain 4: Iteration: 9000 / 10000 [ 90%]  (Sampling)
Chain 1: Iteration: 10000 / 10000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 17.9894 seconds (Warm-up)
Chain 1:                16.7475 seconds (Sampling)
Chain 1:                34.7369 seconds (Total)
Chain 1: 
Chain 2: Iteration: 10000 / 10000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 17.1323 seconds (Warm-up)
Chain 2:                17.7856 seconds (Sampling)
Chain 2:                34.9179 seconds (Total)
Chain 2: 
Chain 3: Iteration: 10000 / 10000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 17.7312 seconds (Warm-up)
Chain 3:                17.3252 seconds (Sampling)
Chain 3:                35.0564 seconds (Total)
Chain 3: 
Chain 4: Iteration: 10000 / 10000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 17.6394 seconds (Warm-up)
Chain 4:                17.244 seconds (Sampling)
Chain 4:                34.8834 seconds (Total)
Chain 4: 

Priors & chain diagnostics

prior_summary(m1d)
Priors for model 'm1d' 
------
Intercept (after predictors centered)
 ~ normal(location = 0, scale = 2.5)

Coefficients
  Specified prior:
    ~ normal(location = [0,0,0,...], scale = [2.5,2.5,2.5,...])
  Adjusted prior:
    ~ normal(location = [0,0,0,...], scale = [5.00,5.00,5.00,...])

Covariance
 ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
------
See help('prior_summary.stanreg') for more details
b1 <- c("(Intercept)", "length_z", "snowt_z", "snowt1_z", "day_z", "bdload_z", "elev_z", "first_c0.353", "male_c0.548", "donor70567", "donor72996")

mcmc_trace(m1d, size = 0.1, pars = b1)

mcmc_dens_overlay(m1d, pars = b1)


summary(m1d)

Model Info:
 function:     stan_glmer
 family:       binomial [logit]
 formula:      surv ~ length_z + snowt_z + snowt1_z + day_z + bdload_z + elev_z + 
       first_c + male_c + donor + (1 | translocation_id)
 algorithm:    sampling
 sample:       20000 (posterior sample size)
 priors:       see help('prior_summary')
 observations: 767
 groups:       translocation_id (24)

Estimates:
                                                  mean   sd   10%   50%   90%
(Intercept)                                      0.3    0.7 -0.6   0.3   1.2 
length_z                                         0.0    0.2 -0.3   0.0   0.3 
snowt_z                                          0.0    0.6 -0.7   0.0   0.7 
snowt1_z                                        -1.2    0.5 -1.9  -1.2  -0.5 
day_z                                            0.1    0.5 -0.6   0.1   0.8 
bdload_z                                        -0.2    0.2 -0.5  -0.2   0.0 
elev_z                                           1.6    0.6  0.8   1.6   2.4 
first_c0.353                                     0.5    0.5 -0.1   0.5   1.1 
male_c0.548                                      0.3    0.2  0.1   0.3   0.6 
donor70567                                      -1.3    0.9 -2.4  -1.3  -0.2 
donor72996                                      -2.2    0.8 -3.2  -2.2  -1.3 
b[(Intercept) translocation_id:70134_1]         -0.6    0.5 -1.3  -0.6   0.1 
b[(Intercept) translocation_id:70370_1]         -0.6    0.8 -1.6  -0.5   0.4 
b[(Intercept) translocation_id:70370_2]          0.0    0.7 -0.9   0.0   0.9 
b[(Intercept) translocation_id:70413_1]         -0.1    0.7 -1.0  -0.1   0.8 
b[(Intercept) translocation_id:70413_2]          0.4    0.7 -0.5   0.4   1.3 
b[(Intercept) translocation_id:70413_3]         -0.3    0.8 -1.3  -0.3   0.6 
b[(Intercept) translocation_id:70414_1]         -0.9    0.8 -2.0  -0.9   0.1 
b[(Intercept) translocation_id:70449_1]          0.1    0.6 -0.7   0.1   0.9 
b[(Intercept) translocation_id:70449_2]          1.6    0.6  0.9   1.6   2.4 
b[(Intercept) translocation_id:70505_1]          0.1    0.5 -0.6   0.1   0.8 
b[(Intercept) translocation_id:70505_2]          0.1    0.6 -0.7   0.1   0.9 
b[(Intercept) translocation_id:70505_3]          0.2    0.7 -0.7   0.2   1.0 
b[(Intercept) translocation_id:70505_4]         -0.1    0.6 -0.9  -0.1   0.7 
b[(Intercept) translocation_id:70550_1]          0.5    0.5 -0.2   0.5   1.2 
b[(Intercept) translocation_id:70550_2]         -0.1    0.6 -0.9  -0.1   0.6 
b[(Intercept) translocation_id:70556_1]         -0.5    0.7 -1.4  -0.5   0.4 
b[(Intercept) translocation_id:70556_2]         -0.8    0.7 -1.8  -0.8   0.0 
b[(Intercept) translocation_id:70619_1]         -0.5    0.6 -1.2  -0.5   0.2 
b[(Intercept) translocation_id:70628_1]          0.8    0.6  0.0   0.8   1.6 
b[(Intercept) translocation_id:70641_1]          0.2    0.7 -0.7   0.2   1.1 
b[(Intercept) translocation_id:70641_2]         -0.3    0.7 -1.1  -0.2   0.6 
b[(Intercept) translocation_id:70641_3]         -0.7    0.7 -1.7  -0.7   0.2 
b[(Intercept) translocation_id:74976_1]          1.4    0.8  0.4   1.3   2.4 
b[(Intercept) translocation_id:74976_2]          0.0    0.7 -0.8   0.0   0.9 
Sigma[translocation_id:(Intercept),(Intercept)]  0.9    0.5  0.4   0.8   1.6 

Fit Diagnostics:
           mean   sd   10%   50%   90%
mean_PPD 0.3    0.0  0.3   0.3   0.3  

The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).

MCMC diagnostics
                                                mcse Rhat n_eff
(Intercept)                                     0.0  1.0  10818
length_z                                        0.0  1.0  25683
snowt_z                                         0.0  1.0  10328
snowt1_z                                        0.0  1.0  12115
day_z                                           0.0  1.0  10672
bdload_z                                        0.0  1.0  28564
elev_z                                          0.0  1.0  11059
first_c0.353                                    0.0  1.0  11342
male_c0.548                                     0.0  1.0  26534
donor70567                                      0.0  1.0  10846
donor72996                                      0.0  1.0  10753
b[(Intercept) translocation_id:70134_1]         0.0  1.0  14706
b[(Intercept) translocation_id:70370_1]         0.0  1.0  15110
b[(Intercept) translocation_id:70370_2]         0.0  1.0  17141
b[(Intercept) translocation_id:70413_1]         0.0  1.0  13738
b[(Intercept) translocation_id:70413_2]         0.0  1.0  14123
b[(Intercept) translocation_id:70413_3]         0.0  1.0  16112
b[(Intercept) translocation_id:70414_1]         0.0  1.0  17969
b[(Intercept) translocation_id:70449_1]         0.0  1.0  12532
b[(Intercept) translocation_id:70449_2]         0.0  1.0  13621
b[(Intercept) translocation_id:70505_1]         0.0  1.0  14823
b[(Intercept) translocation_id:70505_2]         0.0  1.0  18847
b[(Intercept) translocation_id:70505_3]         0.0  1.0  20672
b[(Intercept) translocation_id:70505_4]         0.0  1.0  17878
b[(Intercept) translocation_id:70550_1]         0.0  1.0  11544
b[(Intercept) translocation_id:70550_2]         0.0  1.0  15234
b[(Intercept) translocation_id:70556_1]         0.0  1.0  14110
b[(Intercept) translocation_id:70556_2]         0.0  1.0  14298
b[(Intercept) translocation_id:70619_1]         0.0  1.0  12428
b[(Intercept) translocation_id:70628_1]         0.0  1.0  11195
b[(Intercept) translocation_id:70641_1]         0.0  1.0  11049
b[(Intercept) translocation_id:70641_2]         0.0  1.0  18120
b[(Intercept) translocation_id:70641_3]         0.0  1.0  16196
b[(Intercept) translocation_id:74976_1]         0.0  1.0  14344
b[(Intercept) translocation_id:74976_2]         0.0  1.0  14465
Sigma[translocation_id:(Intercept),(Intercept)] 0.0  1.0   8036
mean_PPD                                        0.0  1.0  19612
log-posterior                                   0.1  1.0   5727

For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).

Posterior predictive check

pp_check(m1d, plotfun = "stat") +
    xlab("Frog survival rate")

Calculate loo

loo_m1d <- loo(m1d, cores = 4, save_psis = TRUE)
print(loo_m1d)

Computed from 20000 by 767 log-likelihood matrix
------
Monte Carlo SE of elpd_loo is 0.0.

All Pareto k estimates are good (k < 0.5).
See help('pareto-k-diagnostic') for details.

Posterior analysis

tidy(m1d, conf.int = TRUE, conf.level = 0.90)

mcmc_areas(m1d, area_method = "equal height", prob = 0.90, pars = b1) +
  geom_vline(xintercept = 0, color = "red", linetype = "dashed", lwd = 0.5) +  
  ggtitle("Model with standardized predictors & translocation_id as grouping variable",
          "Posterior distributions with medians & 90% uncertainty intervals")


mean(bayes_R2(m1d))
[1] 0.322301
# Frog survival by translocation_id
mcmc_areas_ridges(m1d, regex_pars = "translocation_id") +
  geom_vline(xintercept = 0, color = "red", linetype = "dashed", lwd = 0.5) +  
  ggtitle("Frog survival for each translocation_id")

  • Based on loo, model with translocation_id as grouping variable provides slightly better fit than model with site_id as grouping variable.
  • Based on 90% uncertainty intervals, important predictors are similar to those in model without grouping variable: snowt1_z, elev_z, male_c0.548, donor72996.

Model using standardized predictors and 2 grouping variables: translocation_id nested within site_id)

Model specification

m1e <- stan_glmer(
  surv ~ length_z + snowt_z + snowt1_z + day_z + bdload_z + elev_z + first_c + male_c + donor + (1 | translocation_id) + (1 | site_id),
  data = d1,
  family = binomial,
  chains = 4, 
  iter = 5000*2, 
  cores = 4,
  seed = 84735)
starting worker pid=12308 on localhost:11976 at 13:38:14.420
starting worker pid=12337 on localhost:11976 at 13:38:14.541
starting worker pid=12366 on localhost:11976 at 13:38:14.669
starting worker pid=12395 on localhost:11976 at 13:38:14.798

SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.000122 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.22 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 10000 [  0%]  (Warmup)

SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 0.000117 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 1.17 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 10000 [  0%]  (Warmup)

SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 0.00015 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 1.5 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:    1 / 10000 [  0%]  (Warmup)

SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 0.000148 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 1.48 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:    1 / 10000 [  0%]  (Warmup)
Chain 2: Iteration: 1000 / 10000 [ 10%]  (Warmup)
Chain 1: Iteration: 1000 / 10000 [ 10%]  (Warmup)
Chain 3: Iteration: 1000 / 10000 [ 10%]  (Warmup)
Chain 4: Iteration: 1000 / 10000 [ 10%]  (Warmup)
Chain 1: Iteration: 2000 / 10000 [ 20%]  (Warmup)
Chain 2: Iteration: 2000 / 10000 [ 20%]  (Warmup)
Chain 3: Iteration: 2000 / 10000 [ 20%]  (Warmup)
Chain 4: Iteration: 2000 / 10000 [ 20%]  (Warmup)
Chain 1: Iteration: 3000 / 10000 [ 30%]  (Warmup)
Chain 2: Iteration: 3000 / 10000 [ 30%]  (Warmup)
Chain 3: Iteration: 3000 / 10000 [ 30%]  (Warmup)
Chain 4: Iteration: 3000 / 10000 [ 30%]  (Warmup)
Chain 1: Iteration: 4000 / 10000 [ 40%]  (Warmup)
Chain 2: Iteration: 4000 / 10000 [ 40%]  (Warmup)
Chain 3: Iteration: 4000 / 10000 [ 40%]  (Warmup)
Chain 4: Iteration: 4000 / 10000 [ 40%]  (Warmup)
Chain 1: Iteration: 5000 / 10000 [ 50%]  (Warmup)
Chain 1: Iteration: 5001 / 10000 [ 50%]  (Sampling)
Chain 2: Iteration: 5000 / 10000 [ 50%]  (Warmup)
Chain 2: Iteration: 5001 / 10000 [ 50%]  (Sampling)
Chain 3: Iteration: 5000 / 10000 [ 50%]  (Warmup)
Chain 3: Iteration: 5001 / 10000 [ 50%]  (Sampling)
Chain 4: Iteration: 5000 / 10000 [ 50%]  (Warmup)
Chain 4: Iteration: 5001 / 10000 [ 50%]  (Sampling)
Chain 2: Iteration: 6000 / 10000 [ 60%]  (Sampling)
Chain 3: Iteration: 6000 / 10000 [ 60%]  (Sampling)
Chain 1: Iteration: 6000 / 10000 [ 60%]  (Sampling)
Chain 4: Iteration: 6000 / 10000 [ 60%]  (Sampling)
Chain 2: Iteration: 7000 / 10000 [ 70%]  (Sampling)
Chain 3: Iteration: 7000 / 10000 [ 70%]  (Sampling)
Chain 1: Iteration: 7000 / 10000 [ 70%]  (Sampling)
Chain 2: Iteration: 8000 / 10000 [ 80%]  (Sampling)
Chain 3: Iteration: 8000 / 10000 [ 80%]  (Sampling)
Chain 4: Iteration: 7000 / 10000 [ 70%]  (Sampling)
Chain 2: Iteration: 9000 / 10000 [ 90%]  (Sampling)
Chain 3: Iteration: 9000 / 10000 [ 90%]  (Sampling)
Chain 1: Iteration: 8000 / 10000 [ 80%]  (Sampling)
Chain 4: Iteration: 8000 / 10000 [ 80%]  (Sampling)
Chain 2: Iteration: 10000 / 10000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 33.7892 seconds (Warm-up)
Chain 2:                29.5968 seconds (Sampling)
Chain 2:                63.386 seconds (Total)
Chain 2: 
Chain 3: Iteration: 10000 / 10000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 33.485 seconds (Warm-up)
Chain 3:                29.5418 seconds (Sampling)
Chain 3:                63.0268 seconds (Total)
Chain 3: 
Chain 1: Iteration: 9000 / 10000 [ 90%]  (Sampling)
Chain 4: Iteration: 9000 / 10000 [ 90%]  (Sampling)
Chain 1: Iteration: 10000 / 10000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 32.3165 seconds (Warm-up)
Chain 1:                44.1083 seconds (Sampling)
Chain 1:                76.4248 seconds (Total)
Chain 1: 
Chain 4: Iteration: 10000 / 10000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 35.4405 seconds (Warm-up)
Chain 4:                44.1243 seconds (Sampling)
Chain 4:                79.5647 seconds (Total)
Chain 4: 

Priors & chain diagnostics

prior_summary(m1e)
Priors for model 'm1e' 
------
Intercept (after predictors centered)
 ~ normal(location = 0, scale = 2.5)

Coefficients
  Specified prior:
    ~ normal(location = [0,0,0,...], scale = [2.5,2.5,2.5,...])
  Adjusted prior:
    ~ normal(location = [0,0,0,...], scale = [5.00,5.00,5.00,...])

Covariance
 ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
------
See help('prior_summary.stanreg') for more details
b1 <- c("(Intercept)", "length_z", "snowt_z", "snowt1_z", "day_z", "bdload_z", "elev_z", "first_c0.353", "male_c0.548", "donor70567", "donor72996")

mcmc_trace(m1e, size = 0.1, pars = b1)

mcmc_dens_overlay(m1e, pars = b1)


summary(m1e)

Model Info:
 function:     stan_glmer
 family:       binomial [logit]
 formula:      surv ~ length_z + snowt_z + snowt1_z + day_z + bdload_z + elev_z + 
       first_c + male_c + donor + (1 | translocation_id) + (1 | 
       site_id)
 algorithm:    sampling
 sample:       20000 (posterior sample size)
 priors:       see help('prior_summary')
 observations: 767
 groups:       translocation_id (24), site_id (12)

Estimates:
                                                  mean   sd   10%   50%   90%
(Intercept)                                      0.1    1.0 -1.2   0.1   1.3 
length_z                                         0.0    0.2 -0.3   0.0   0.3 
snowt_z                                         -0.1    0.5 -0.7  -0.1   0.6 
snowt1_z                                        -0.5    0.7 -1.5  -0.6   0.4 
day_z                                            0.2    0.5 -0.5   0.2   0.9 
bdload_z                                        -0.2    0.2 -0.5  -0.2   0.0 
elev_z                                           1.7    0.9  0.6   1.7   2.9 
first_c0.353                                     0.4    0.4  0.0   0.4   0.9 
male_c0.548                                      0.3    0.2  0.1   0.3   0.6 
donor70567                                      -1.0    1.5 -2.7  -1.1   0.8 
donor72996                                      -2.0    1.1 -3.4  -2.1  -0.7 
b[(Intercept) translocation_id:70134_1]         -0.2    0.5 -0.9  -0.1   0.4 
b[(Intercept) translocation_id:70370_1]         -0.4    0.6 -1.1  -0.3   0.2 
b[(Intercept) translocation_id:70370_2]          0.1    0.5 -0.5   0.1   0.7 
b[(Intercept) translocation_id:70413_1]          0.0    0.5 -0.6   0.0   0.6 
b[(Intercept) translocation_id:70413_2]          0.2    0.5 -0.4   0.1   0.9 
b[(Intercept) translocation_id:70413_3]         -0.2    0.6 -0.9  -0.2   0.4 
b[(Intercept) translocation_id:70414_1]         -0.3    0.7 -1.2  -0.2   0.3 
b[(Intercept) translocation_id:70449_1]         -0.1    0.5 -0.7  -0.1   0.5 
b[(Intercept) translocation_id:70449_2]          0.7    0.7  0.0   0.6   1.7 
b[(Intercept) translocation_id:70505_1]          0.0    0.5 -0.5   0.0   0.6 
b[(Intercept) translocation_id:70505_2]          0.1    0.5 -0.5   0.1   0.7 
b[(Intercept) translocation_id:70505_3]          0.1    0.5 -0.5   0.0   0.7 
b[(Intercept) translocation_id:70505_4]         -0.1    0.5 -0.7  -0.1   0.5 
b[(Intercept) translocation_id:70550_1]          0.3    0.5 -0.3   0.2   0.9 
b[(Intercept) translocation_id:70550_2]         -0.1    0.5 -0.7  -0.1   0.4 
b[(Intercept) translocation_id:70556_1]         -0.2    0.6 -0.9  -0.1   0.4 
b[(Intercept) translocation_id:70556_2]         -0.2    0.6 -1.0  -0.1   0.4 
b[(Intercept) translocation_id:70619_1]         -0.2    0.5 -0.9  -0.1   0.4 
b[(Intercept) translocation_id:70628_1]          0.3    0.6 -0.3   0.2   1.1 
b[(Intercept) translocation_id:70641_1]          0.2    0.5 -0.4   0.1   0.9 
b[(Intercept) translocation_id:70641_2]         -0.1    0.5 -0.7   0.0   0.5 
b[(Intercept) translocation_id:70641_3]         -0.3    0.6 -1.1  -0.2   0.2 
b[(Intercept) translocation_id:74976_1]          0.5    0.7 -0.1   0.4   1.5 
b[(Intercept) translocation_id:74976_2]         -0.1    0.5 -0.7  -0.1   0.5 
b[(Intercept) site_id:70134]                    -0.3    0.6 -1.1  -0.3   0.4 
b[(Intercept) site_id:70370]                    -0.9    1.0 -2.3  -0.7   0.2 
b[(Intercept) site_id:70413]                     0.0    1.1 -1.3   0.0   1.2 
b[(Intercept) site_id:70414]                    -0.7    1.1 -2.1  -0.5   0.3 
b[(Intercept) site_id:70449]                     1.0    0.8  0.0   1.0   2.0 
b[(Intercept) site_id:70505]                     0.2    0.6 -0.5   0.1   0.9 
b[(Intercept) site_id:70550]                     0.3    0.6 -0.4   0.2   1.1 
b[(Intercept) site_id:70556]                    -0.7    0.9 -1.8  -0.5   0.3 
b[(Intercept) site_id:70619]                    -0.4    0.7 -1.3  -0.3   0.4 
b[(Intercept) site_id:70628]                     0.6    0.8 -0.2   0.5   1.7 
b[(Intercept) site_id:70641]                    -0.1    0.7 -1.0  -0.1   0.7 
b[(Intercept) site_id:74976]                     0.7    0.9 -0.2   0.6   1.9 
Sigma[translocation_id:(Intercept),(Intercept)]  0.4    0.4  0.0   0.3   0.9 
Sigma[site_id:(Intercept),(Intercept)]           1.2    1.3  0.1   0.8   2.6 

Fit Diagnostics:
           mean   sd   10%   50%   90%
mean_PPD 0.3    0.0  0.3   0.3   0.3  

The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).

MCMC diagnostics
                                                mcse Rhat n_eff
(Intercept)                                     0.0  1.0   8760
length_z                                        0.0  1.0  26927
snowt_z                                         0.0  1.0  12055
snowt1_z                                        0.0  1.0   5542
day_z                                           0.0  1.0  10670
bdload_z                                        0.0  1.0  28589
elev_z                                          0.0  1.0  10334
first_c0.353                                    0.0  1.0  12870
male_c0.548                                     0.0  1.0  27066
donor70567                                      0.0  1.0   9522
donor72996                                      0.0  1.0   8564
b[(Intercept) translocation_id:70134_1]         0.0  1.0  10178
b[(Intercept) translocation_id:70370_1]         0.0  1.0  13784
b[(Intercept) translocation_id:70370_2]         0.0  1.0  21315
b[(Intercept) translocation_id:70413_1]         0.0  1.0  18159
b[(Intercept) translocation_id:70413_2]         0.0  1.0  16405
b[(Intercept) translocation_id:70413_3]         0.0  1.0  18489
b[(Intercept) translocation_id:70414_1]         0.0  1.0   8310
b[(Intercept) translocation_id:70449_1]         0.0  1.0  15651
b[(Intercept) translocation_id:70449_2]         0.0  1.0   4203
b[(Intercept) translocation_id:70505_1]         0.0  1.0  18690
b[(Intercept) translocation_id:70505_2]         0.0  1.0  23587
b[(Intercept) translocation_id:70505_3]         0.0  1.0  23491
b[(Intercept) translocation_id:70505_4]         0.0  1.0  22403
b[(Intercept) translocation_id:70550_1]         0.0  1.0  12418
b[(Intercept) translocation_id:70550_2]         0.0  1.0  19758
b[(Intercept) translocation_id:70556_1]         0.0  1.0  14624
b[(Intercept) translocation_id:70556_2]         0.0  1.0   7455
b[(Intercept) translocation_id:70619_1]         0.0  1.0  11194
b[(Intercept) translocation_id:70628_1]         0.0  1.0   8502
b[(Intercept) translocation_id:70641_1]         0.0  1.0  15452
b[(Intercept) translocation_id:70641_2]         0.0  1.0  20501
b[(Intercept) translocation_id:70641_3]         0.0  1.0   9364
b[(Intercept) translocation_id:74976_1]         0.0  1.0   5537
b[(Intercept) translocation_id:74976_2]         0.0  1.0  20547
b[(Intercept) site_id:70134]                    0.0  1.0  12998
b[(Intercept) site_id:70370]                    0.0  1.0   5422
b[(Intercept) site_id:70413]                    0.0  1.0  14150
b[(Intercept) site_id:70414]                    0.0  1.0   9926
b[(Intercept) site_id:70449]                    0.0  1.0   4422
b[(Intercept) site_id:70505]                    0.0  1.0  10980
b[(Intercept) site_id:70550]                    0.0  1.0   9459
b[(Intercept) site_id:70556]                    0.0  1.0  10122
b[(Intercept) site_id:70619]                    0.0  1.0   9424
b[(Intercept) site_id:70628]                    0.0  1.0   6708
b[(Intercept) site_id:70641]                    0.0  1.0   9372
b[(Intercept) site_id:74976]                    0.0  1.0   7086
Sigma[translocation_id:(Intercept),(Intercept)] 0.0  1.0   4082
Sigma[site_id:(Intercept),(Intercept)]          0.0  1.0   4867
mean_PPD                                        0.0  1.0  19556
log-posterior                                   0.1  1.0   5917

For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).

Posterior predictive check

pp_check(m1e, plotfun = "stat") +
    xlab("Frog survival rate")

Calculate loo

loo_m1e <- loo(m1e, cores = 4, save_psis = TRUE)
print(loo_m1e)

Computed from 20000 by 767 log-likelihood matrix
------
Monte Carlo SE of elpd_loo is 0.0.

All Pareto k estimates are good (k < 0.5).
See help('pareto-k-diagnostic') for details.

Posterior analysis

tidy(m1e, conf.int = TRUE, conf.level = 0.90)

mcmc_areas(m1e, area_method = "equal height", prob = 0.90, pars = b1) +
  geom_vline(xintercept = 0, color = "red", linetype = "dashed", lwd = 0.5) +  
  ggtitle("Model with standardized predictors & translocation_id as grouping variable",
          "Posterior distributions with medians & 90% uncertainty intervals")


mean(bayes_R2(m1e))
[1] 0.3232688
# Frog survival by translocation_id
mcmc_areas_ridges(m1e, regex_pars = ("translocation_id")) +
  geom_vline(xintercept = 0, color = "red", linetype = "dashed", lwd = 0.5) +  
  ggtitle("Frog survival for each translocation_id")


# Frog survival by site_id
mcmc_areas_ridges(m1e, regex_pars = ("site_id")) +
  geom_vline(xintercept = 0, color = "red", linetype = "dashed", lwd = 0.5) +  
  ggtitle("Frog survival for each translocation_id")

Compare 4 models with different grouping structures

loo_compare(loo_m1b, loo_m1c, loo_m1d, loo_m1e)
    elpd_diff se_diff
m1c   0.0       0.0  
m1e  -0.5       1.9  
m1d  -1.1       3.1  
m1b -18.2       6.3  

Additional analysis of model m1d

Posterior classification accuracy

# Posterior predictive models for each frog in dataset
set.seed(84735)
survival_predict_1 <- posterior_predict(m1d, newdata = d1)

# Classify frog survival based on mean predicted survival probability
survival_predict_2 <- d1 %>% 
  mutate(survival_prob = colMeans(survival_predict_1),
         survival_class = as.numeric(survival_prob >= 0.5)) %>% 
  select(site_id, translocation_id, year, snowt1_z, elev_z, male_c, donor, surv, survival_prob, survival_class)

# Generate confusion matrix and calculate accuracy
survival_predict_2 %>% 
  tabyl(surv, survival_class) %>% 
  adorn_totals(c("row", "col"))
  surv   0   1 Total
     0 471  62   533
     1  96 138   234
 Total 567 200   767
(471 + 138)/767 # overall_accuracy 
[1] 0.7940026
471/533 # true negative rate
[1] 0.8836773
138/234 # true positive rate
[1] 0.5897436
# Generate confusion matrix using function
set.seed(84735)
classification_summary(m1d, data = d1, cutoff = 0.5)
$confusion_matrix
 y   0   1
 0 471  62
 1  96 138

$accuracy_rates
NA
  • Sensitivity or “true positive rate” measures the proportion of Y = 1 observations that are accurately classified.
  • Specificity or “true negative rate” measures the proportion of Y = 0 observations that are accurately classified (i.e., proportion of frogs correctly identified as not surviving).

Create reduced model containing only important predictors for subsequent analyses

  • Reduced model is helpful because it reduces the number of values used to calculate predictive models.
m1d_reduce <- stan_glmer(
  surv ~ snowt1_z + elev_z + male_c + donor + (1 | translocation_id),
  data = d1,
  family = binomial,
  chains = 4, 
  iter = 5000*2, 
  cores = 4,
  seed = 84735)
starting worker pid=12738 on localhost:11976 at 13:40:15.905
starting worker pid=12767 on localhost:11976 at 13:40:16.051
starting worker pid=12796 on localhost:11976 at 13:40:16.183
starting worker pid=12825 on localhost:11976 at 13:40:16.313

SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.000102 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.02 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 10000 [  0%]  (Warmup)

SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 0.000101 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 1.01 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 10000 [  0%]  (Warmup)

SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 0.000108 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 1.08 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:    1 / 10000 [  0%]  (Warmup)

SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 0.000139 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 1.39 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:    1 / 10000 [  0%]  (Warmup)
Chain 1: Iteration: 1000 / 10000 [ 10%]  (Warmup)
Chain 3: Iteration: 1000 / 10000 [ 10%]  (Warmup)
Chain 2: Iteration: 1000 / 10000 [ 10%]  (Warmup)
Chain 4: Iteration: 1000 / 10000 [ 10%]  (Warmup)
Chain 1: Iteration: 2000 / 10000 [ 20%]  (Warmup)
Chain 2: Iteration: 2000 / 10000 [ 20%]  (Warmup)
Chain 3: Iteration: 2000 / 10000 [ 20%]  (Warmup)
Chain 4: Iteration: 2000 / 10000 [ 20%]  (Warmup)
Chain 1: Iteration: 3000 / 10000 [ 30%]  (Warmup)
Chain 2: Iteration: 3000 / 10000 [ 30%]  (Warmup)
Chain 3: Iteration: 3000 / 10000 [ 30%]  (Warmup)
Chain 4: Iteration: 3000 / 10000 [ 30%]  (Warmup)
Chain 1: Iteration: 4000 / 10000 [ 40%]  (Warmup)
Chain 2: Iteration: 4000 / 10000 [ 40%]  (Warmup)
Chain 3: Iteration: 4000 / 10000 [ 40%]  (Warmup)
Chain 4: Iteration: 4000 / 10000 [ 40%]  (Warmup)
Chain 1: Iteration: 5000 / 10000 [ 50%]  (Warmup)
Chain 1: Iteration: 5001 / 10000 [ 50%]  (Sampling)
Chain 2: Iteration: 5000 / 10000 [ 50%]  (Warmup)
Chain 2: Iteration: 5001 / 10000 [ 50%]  (Sampling)
Chain 3: Iteration: 5000 / 10000 [ 50%]  (Warmup)
Chain 3: Iteration: 5001 / 10000 [ 50%]  (Sampling)
Chain 4: Iteration: 5000 / 10000 [ 50%]  (Warmup)
Chain 4: Iteration: 5001 / 10000 [ 50%]  (Sampling)
Chain 1: Iteration: 6000 / 10000 [ 60%]  (Sampling)
Chain 2: Iteration: 6000 / 10000 [ 60%]  (Sampling)
Chain 3: Iteration: 6000 / 10000 [ 60%]  (Sampling)
Chain 4: Iteration: 6000 / 10000 [ 60%]  (Sampling)
Chain 1: Iteration: 7000 / 10000 [ 70%]  (Sampling)
Chain 2: Iteration: 7000 / 10000 [ 70%]  (Sampling)
Chain 3: Iteration: 7000 / 10000 [ 70%]  (Sampling)
Chain 4: Iteration: 7000 / 10000 [ 70%]  (Sampling)
Chain 1: Iteration: 8000 / 10000 [ 80%]  (Sampling)
Chain 2: Iteration: 8000 / 10000 [ 80%]  (Sampling)
Chain 3: Iteration: 8000 / 10000 [ 80%]  (Sampling)
Chain 4: Iteration: 8000 / 10000 [ 80%]  (Sampling)
Chain 2: Iteration: 9000 / 10000 [ 90%]  (Sampling)
Chain 1: Iteration: 9000 / 10000 [ 90%]  (Sampling)
Chain 3: Iteration: 9000 / 10000 [ 90%]  (Sampling)
Chain 4: Iteration: 9000 / 10000 [ 90%]  (Sampling)
Chain 2: Iteration: 10000 / 10000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 12.8679 seconds (Warm-up)
Chain 2:                11.9961 seconds (Sampling)
Chain 2:                24.864 seconds (Total)
Chain 2: 
Chain 1: Iteration: 10000 / 10000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 12.7424 seconds (Warm-up)
Chain 1:                12.4851 seconds (Sampling)
Chain 1:                25.2275 seconds (Total)
Chain 1: 
Chain 3: Iteration: 10000 / 10000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 12.7796 seconds (Warm-up)
Chain 3:                12.5353 seconds (Sampling)
Chain 3:                25.3149 seconds (Total)
Chain 3: 
Chain 4: Iteration: 10000 / 10000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 13.0358 seconds (Warm-up)
Chain 4:                12.7509 seconds (Sampling)
Chain 4:                25.7867 seconds (Total)
Chain 4: 

Confirm fit of reduced model is comparable to full model

loo_m1d_reduce <- loo(m1d, cores = 4, save_psis = TRUE)
print(loo_m1d)

Computed from 20000 by 767 log-likelihood matrix
------
Monte Carlo SE of elpd_loo is 0.0.

All Pareto k estimates are good (k < 0.5).
See help('pareto-k-diagnostic') for details.
loo_compare(loo_m1d, loo_m1d_reduce)
    elpd_diff se_diff
m1d 0.0       0.0    
m1d 0.0       0.0    

Simulate and plot posterior predictive models & observations

m1d_group_means <- d1 %>% 
  group_by(translocation_id) %>% 
  summarize(count = n(), psurv = mean(surv), elev_z = mean(elev_z), snowt1_z = mean(snowt1_z)) %>% 
  mutate(
    male_c = as.factor(0.548), 
    donor = case_when(
      str_detect(translocation_id, "74976") | str_detect(translocation_id, "70556") ~ "70459",
      str_detect(translocation_id, "70413") ~ "70567",
      TRUE ~ "72996")) %>% 
  arrange(psurv)  # to arrange plot values/labels by psurv

predictions_m1d <- posterior_predict(m1d_reduce, newdata = m1d_group_means)

ppc_intervals(
  y = m1d_group_means$psurv,
  yrep = predictions_m1d, 
  prob = 0.5, 
  prob_outer = 0.8) +
ggplot2::scale_x_continuous(
  labels = m1d_group_means$translocation_id,
  breaks = 1:nrow(m1d_group_means)) +
  xaxis_text(angle = 90, vjust = 0.5) +
xlab("Translocation_id")

  • Uncertainty intervals appear incorrect for binomial model, but not sure how to adjust.
  • Based on mean values, model does an adequate job of predicting survival based on group-level and individual-level characteristics.
  • What explains the low predicted survival at 70370 (both translocations), and high versus low predicted survival at 70550 for 1st and 2nd translocations, respectively? Both sites seem to contain high quality habitat (i.e., high elevation).
  • The predicted survival is based on assumption that only males are included. Run another set of models for females?

Plot posterior plausible relationships between each predictor and the probability of frog survival (conditional effects)

  • As described in https://rdrr.io/cran/brms/man/conditional_effects.brmsfit.html, “When creating conditional_effects for a particular predictor (or interaction of two predictors), one has to choose the values of all other predictors to condition on. By default, the mean is used for continuous variables and the reference category is used for factors, but you may change these values via argument conditions. This also has an implication for the points argument: In the created plots, only those points will be shown that correspond to the factor levels actually used in the conditioning, in order not to create the false impression of bad model fit, where it is just due to conditioning on certain factor levels.”
  • Using rstanarm, I was unable to do this as I expected to: Start with data = d1 and using “newdata =” to specify that predictions should be based on data in d1_elevz (for example), in which non-plotted continuous predictors are set to the mean value. Not sure if my alternative approach is correct.
  • Approach used: Instead of starting the routine with d1, I started with the new dataset in which all variables except the one of interest were set to mean values (continuous variables) and the reference level (categorical variables). This produced plots/results that were very similar to those produced in brms using the conditional_effects function.

Effect of elevation

d1_elevz <- d1 %>% 
  mutate(snowt1_z = mean(snowt1_z), 
         male_c = recode(male_c, "-0.452" = "0.548"),
         donor = recode(donor, "70567" = "70459", "72996" = "70459"))

d1_elevz %>% distinct(snowt1_z, male_c, donor) # ensure variables were set to correct values
  
d1_elevz %>% 
  select(site_id, translocation_id, snowt1_z, elev_z, male_c, donor, surv) %>% 
  add_epred_draws(m1d_reduce, ndraws = 100, re_formula = NA) %>% 
  ggplot(aes(x = elev_z, y = surv)) +
    geom_line(aes(y = .epred, group = .draw), color = "blue", size = 0.2, alpha = 0.5) +
#    geom_point(aes(y = .epred)) +
  labs(x = "Elevation of recipient site (standardized)", y = "Probability of 1-year frog survival")



# d1 %>%  # Test of add_predicted_draws used for similar purpose - unsuccessful
#   select(site_id, translocation_id, snowt1_z, elev_z, male_c, donor, surv) %>% 
#   add_predicted_draws(m1d_reduce, ndraws = 100, re_formula = NA, newdata = d1_elevz) %>% View()
#   ggplot(aes(x = elev_z, y = surv)) +
#     geom_line(aes(y = .epred, group = .draw), color = "blue", size = 0.2, alpha = 0.5) +
#     geom_point(aes(y = .epred)) +
#   labs(x = "Elevation of recipient site (standardized)", y = "Probability of 1-year frog survival")

Effect of winter severity in year t1

d1_snowt1z <- d1 %>% 
  mutate(elev_z = mean(elev_z), 
         male_c = recode(male_c, "-0.452" = "0.548"),
         donor = recode(donor, "70567" = "70459", "72996" = "70459"))

d1_snowt1z %>% distinct(elev_z, male_c, donor) # ensure variables were set to correct values
  
d1_snowt1z %>% 
  select(site_id, translocation_id, snowt1_z, elev_z, male_c, donor, surv) %>% 
  add_epred_draws(m1d_reduce, ndraws = 100, re_formula = NA) %>% 
  ggplot(aes(x = snowt1_z, y = surv)) +
    geom_line(aes(y = .epred, group = .draw), color = "blue", size = 0.2, alpha = 0.5) +
#    geom_point(aes(y = .epred)) +
  labs(x = "Winter severity in year following translocation (standardized)", y = "Probability of 1-year frog survival") 

Effect of sex

d1_malec <- d1 %>% 
  mutate(elev_z = mean(elev_z), 
         snowt1_z = mean(snowt1_z),
         donor = recode(donor, "70567" = "70459", "72996" = "70459"))

d1_malec %>% distinct(elev_z, snowt1_z, donor) # ensure variables were set to correct values
  
d1_malec %>% 
  select(site_id, translocation_id, snowt1_z, elev_z, male_c, donor, surv) %>% 
  add_epred_draws(m1d_reduce, ndraws = 100, re_formula = NA) %>% 
  ggplot(aes(x = male_c, y = surv)) +
    geom_boxplot(aes(y = .epred, x = male_c)) +
#    geom_point(aes(y = .epred)) +
    labs(x = "Sex (centered)", y = "Probability of 1-year frog survival")

Effect of donor

d1_donor <- d1 %>% 
  mutate(elev_z = mean(elev_z), 
         snowt1_z = mean(snowt1_z),
         male_c = recode(male_c, "-0.452" = "0.548"))

d1_donor %>% distinct(elev_z, snowt1_z, male_c) # ensure variables were set to correct values
  
d1_donor %>% 
  select(site_id, translocation_id, snowt1_z, elev_z, male_c, donor, surv) %>% 
  add_epred_draws(m1d_reduce, ndraws = 100, re_formula = NA) %>% 
  ggplot(aes(x = donor, y = surv)) +
    geom_boxplot(aes(y = .epred, x = donor)) +
#    geom_point(aes(y = .epred)) +
    labs(x = "Donor population", y = "Probability of 1-year frog survival")

Joint effects of elevation and sex

d1_elevz_malec <- d1 %>% 
  mutate(snowt1_z = mean(snowt1_z), 
         donor = recode(donor, "70567" = "70459", "72996" = "70459"))

d1_elevz_malec %>% distinct(snowt1_z, donor) # ensure variables were set to correct values
  
pal <- c("#0d0887", "#f89540")
d1_elevz_malec %>% 
  select(site_id, translocation_id, snowt1_z, elev_z, male_c, donor, surv) %>% 
  add_epred_draws(m1d_reduce, ndraws = 100, re_formula = NA) %>% 
  ggplot(aes(x = elev_z, y = surv, color = male_c)) +
    geom_line(aes(y = .epred, group = paste(male_c, .draw)), size = 0.4, alpha = 0.5) +
#    geom_point(aes(y = .epred)) +
    labs(x = "Elevation of recipient site (standardized)", y = "Probability of 1-year frog survival") +
    scale_color_manual(values = pal)  + 
    theme_classic()

For comparison, fit model using brms and conditional_effects function to plot variable-specific effects

Fit model

m1d_reduce_brms <- brm(
  surv ~ snowt1_z + elev_z + male_c + donor + (1 | translocation_id),
  data = d1,
  family = bernoulli(),
  chains = 4, 
  iter = 5000*2, 
  cores = 4,
  seed = 84735)
Compiling Stan program...
Start sampling
starting worker pid=13412 on localhost:11976 at 13:41:34.174
starting worker pid=13441 on localhost:11976 at 13:41:34.302
starting worker pid=13470 on localhost:11976 at 13:41:34.439
starting worker pid=13499 on localhost:11976 at 13:41:34.568

SAMPLING FOR MODEL 'dbb8efe82731517ccb01114eea2991f6' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 4.4e-05 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.44 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 10000 [  0%]  (Warmup)
Chain 1: Iteration: 1000 / 10000 [ 10%]  (Warmup)

SAMPLING FOR MODEL 'dbb8efe82731517ccb01114eea2991f6' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 4.7e-05 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.47 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 10000 [  0%]  (Warmup)
Chain 1: Iteration: 2000 / 10000 [ 20%]  (Warmup)
Chain 2: Iteration: 1000 / 10000 [ 10%]  (Warmup)
Chain 1: Iteration: 3000 / 10000 [ 30%]  (Warmup)

SAMPLING FOR MODEL 'dbb8efe82731517ccb01114eea2991f6' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 7.9e-05 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.79 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:    1 / 10000 [  0%]  (Warmup)
Chain 2: Iteration: 2000 / 10000 [ 20%]  (Warmup)
Chain 1: Iteration: 4000 / 10000 [ 40%]  (Warmup)

SAMPLING FOR MODEL 'dbb8efe82731517ccb01114eea2991f6' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 4.9e-05 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.49 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:    1 / 10000 [  0%]  (Warmup)
Chain 2: Iteration: 3000 / 10000 [ 30%]  (Warmup)
Chain 1: Iteration: 5000 / 10000 [ 50%]  (Warmup)
Chain 1: Iteration: 5001 / 10000 [ 50%]  (Sampling)
Chain 3: Iteration: 1000 / 10000 [ 10%]  (Warmup)
Chain 4: Iteration: 1000 / 10000 [ 10%]  (Warmup)
Chain 2: Iteration: 4000 / 10000 [ 40%]  (Warmup)
Chain 3: Iteration: 2000 / 10000 [ 20%]  (Warmup)
Chain 1: Iteration: 6000 / 10000 [ 60%]  (Sampling)
Chain 4: Iteration: 2000 / 10000 [ 20%]  (Warmup)
Chain 2: Iteration: 5000 / 10000 [ 50%]  (Warmup)
Chain 2: Iteration: 5001 / 10000 [ 50%]  (Sampling)
Chain 3: Iteration: 3000 / 10000 [ 30%]  (Warmup)
Chain 4: Iteration: 3000 / 10000 [ 30%]  (Warmup)
Chain 1: Iteration: 7000 / 10000 [ 70%]  (Sampling)
Chain 3: Iteration: 4000 / 10000 [ 40%]  (Warmup)
Chain 4: Iteration: 4000 / 10000 [ 40%]  (Warmup)
Chain 2: Iteration: 6000 / 10000 [ 60%]  (Sampling)
Chain 3: Iteration: 5000 / 10000 [ 50%]  (Warmup)
Chain 3: Iteration: 5001 / 10000 [ 50%]  (Sampling)
Chain 1: Iteration: 8000 / 10000 [ 80%]  (Sampling)
Chain 4: Iteration: 5000 / 10000 [ 50%]  (Warmup)
Chain 4: Iteration: 5001 / 10000 [ 50%]  (Sampling)
Chain 2: Iteration: 7000 / 10000 [ 70%]  (Sampling)
Chain 3: Iteration: 6000 / 10000 [ 60%]  (Sampling)
Chain 1: Iteration: 9000 / 10000 [ 90%]  (Sampling)
Chain 4: Iteration: 6000 / 10000 [ 60%]  (Sampling)
Chain 3: Iteration: 7000 / 10000 [ 70%]  (Sampling)
Chain 2: Iteration: 8000 / 10000 [ 80%]  (Sampling)
Chain 1: Iteration: 10000 / 10000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 3.11244 seconds (Warm-up)
Chain 1:                3.90152 seconds (Sampling)
Chain 1:                7.01396 seconds (Total)
Chain 1: 
Chain 4: Iteration: 7000 / 10000 [ 70%]  (Sampling)
Chain 3: Iteration: 8000 / 10000 [ 80%]  (Sampling)
Chain 2: Iteration: 9000 / 10000 [ 90%]  (Sampling)
Chain 3: Iteration: 9000 / 10000 [ 90%]  (Sampling)
Chain 4: Iteration: 8000 / 10000 [ 80%]  (Sampling)
Chain 2: Iteration: 10000 / 10000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 3.07934 seconds (Warm-up)
Chain 2:                4.26893 seconds (Sampling)
Chain 2:                7.34827 seconds (Total)
Chain 2: 
Chain 3: Iteration: 10000 / 10000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 3.08411 seconds (Warm-up)
Chain 3:                3.02089 seconds (Sampling)
Chain 3:                6.105 seconds (Total)
Chain 3: 
Chain 4: Iteration: 9000 / 10000 [ 90%]  (Sampling)
Chain 4: Iteration: 10000 / 10000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 2.96916 seconds (Warm-up)
Chain 4:                4.03496 seconds (Sampling)
Chain 4:                7.00411 seconds (Total)
Chain 4: 

Plot conditional effects

plot(conditional_effects(m1d_reduce_brms, spaghetti = TRUE, mean = TRUE, ndraws = 100), ask = FALSE)

TO DO

LS0tCnRpdGxlOiAiUHJlZGljdG9ycyBvZiBmcm9nIHN1cnZpdmFsIGZvbGxvd2luZyB0cmFuc2xvY2F0aW9ucyAtIGFuYWx5c2VzIgphdXRob3I6ICJSb2xhbmQgS25hcHAiCm91dHB1dDogaHRtbF9ub3RlYm9vawotLS0KClRoZSBnb2FsIG9mIHRoZSBhbmFseXNlcyBpbiB0aGlzIG5vdGVib29rIGlzIHRvIGlkZW50aWZ5IHByZWRpY3RvcnMgb2YgZnJvZyBzdXJ2aXZhbCBmb2xsb3dpbmcgdHJhbnNsb2NhdGlvbi4gVGhlIGFuYWx5c2VzIGFyZSBjb25kdWN0ZWQgaW4gYSBCYXllc2lhbiBmcmFtZXdvcmsgdXNpbmcgdGhlIFIgcGFja2FnZSByc3RhbmFybS4KCiMjIExvYWQgcGFja2FnZXMKCmBgYHtyIGxvYWQtcGFja2FnZXN9CmxpYnJhcnkodGlkeXZlcnNlKQpsaWJyYXJ5KGJybXMpCmxpYnJhcnkocnN0YW5hcm0pCmxpYnJhcnkoYmF5ZXNwbG90KQpsaWJyYXJ5KHRpZHliYXllcykKbGlicmFyeShsb28pCmxpYnJhcnkoYnJvb20ubWl4ZWQpCmxpYnJhcnkoamFuaXRvcikKCnNvdXJjZSgiY2xhc3NpZmljYXRpb25fc3VtbWFyeS5SIikKYGBgCgojIyBSZWFkIGluIGZyb2cgdHJhbnNsb2NhdGlvbiBkYXRhc2V0ICAKKiBSZWFkIGluIGRhdGEsIGNyZWF0ZSB2YXJpYWJsZXMgZm9yIGFuYWx5c2lzLCBzdGFuZGFyZGl6ZSBjb250aW51b3VzIHZhcmlhYmxlcy4gCiogYXJtIHBhY2thZ2UgdXNlZCB0byByZXNjYWxlIGNvbnRpbnVvdXMgdmFyaWFibGVzIChzdWJzdHJhY3QgbWVhbiwgZGl2aWRlIGJ5IDIgU0QpLCBtYW51YWwgY2VudGVyaW5nIG9mIGJpbmFyeSB2YXJpYWJsZXMuCmBgYHtyIHJlYWQtdHJhbnNsb2NhdGlvbi1kYXRhfQpkMSA8LSByZWFkX2NzdihoZXJlOjpoZXJlKCJkYXRhIiwgImNsZWFuIiwgImZyb2dfdHJhbnNsb2NhdGlvbl9maW5hbC5jc3YiKSkgJT4lIAogIGRyb3BfbmEoYmRfbG9hZCkgJT4lICAjIGRyb3AgeHggcmVjb3JkcyBmcm9tIDcweHh4IHdoZXJlIGJkX2xvYWQgaXMgTkEKICBtdXRhdGUoCiAgICBmaXJzdCA9IGlmX2Vsc2Uob3JkZXIgPT0gMSwgMSwgMCksCiAgICBzdXJ2ID0gaWZfZWxzZShzdXJ2aXZhbCA8IDAuNSwgMCwgMSksCiAgICBtYWxlID0gaWZfZWxzZShzZXggPT0gIm0iLCAxLCAwKSwKICAgIGVsZXZfbG8gPSBpZl9lbHNlKGVsZXZhdGlvbiA8IDMwMjAsIDEsIDApLCAjIHVzaW5nIGludGVydmFsIGRlZmluZWQgYnkgY3V0X2ludGVydmFsLCBuID0gMgogICAgZWxldl96ID0gYXJtOjpyZXNjYWxlKGVsZXZhdGlvbiksCiAgICBzbm93dF96ID0gYXJtOjpyZXNjYWxlKHNub3dfdCksCiAgICBzbm93dDFfeiA9IGFybTo6cmVzY2FsZShzbm93X3QxKSwKICAgIGRheV96ID0gYXJtOjpyZXNjYWxlKGRheSksCiAgICBsZW5ndGhfeiA9IGFybTo6cmVzY2FsZShsZW5ndGgpLAogICAgYmRsb2FkX2wgPSBsb2cxMChiZF9sb2FkICsgMSksCiAgICBiZGxvYWRfeiA9IGFybTo6cmVzY2FsZShiZGxvYWRfbCksCiAgICBzaG9yZV9jID0gc2hvcmUgLSBtZWFuKHNob3JlKSwKICAgIG1hbGVfYyA9IG1hbGUgLSBtZWFuKG1hbGUpLAogICAgZWxldmxvX2MgPSBlbGV2X2xvIC0gbWVhbihlbGV2X2xvKSwKICAgIGZpcnN0X2MgPSBmaXJzdCAtIG1lYW4oZmlyc3QpLAogICAgYWNyb3NzKGMoc2hvcmVfYywgbWFsZV9jLCBlbGV2bG9fYywgZmlyc3RfYyksIHJvdW5kLCAzKSwgCiAgICBhY3Jvc3MoYyhzaXRlX2lkLCBzaG9yZV9jLCBmaXJzdF9jLCBkb25vciwgbWFsZV9jLCBlbGV2bG9fYywgeWVhciksIGFzLmZhY3RvciksCiAgICBzdXJ2ID0gYXMuaW50ZWdlcihzdXJ2KSkKCiMgQWRkIGEgZGVzY3JpcHRpdmUgdHJhbnNsb2NhdGlvbl9pZApkMSA8LSBkMSAlPiUgCiAgZGlzdGluY3Qoc2l0ZV9pZCwgeWVhcikgJT4lIAogIGFycmFuZ2Uoc2l0ZV9pZCwgeWVhcikgJT4lIAogIGdyb3VwX2J5KHNpdGVfaWQpICU+JSAKICBtdXRhdGUodHJhbnNubyA9IHNlcV9sZW4obigpKSkgJT4lIAogIHVuZ3JvdXAoc2l0ZV9pZCkgJT4lIAogIHVuaXRlKHRyYW5zbG9jYXRpb25faWQsIGMoInNpdGVfaWQiLCAidHJhbnNubyIpLCByZW1vdmUgPSBGQUxTRSkgJT4lIAogIHNlbGVjdCgtdHJhbnNubykgJT4lIAogIG11dGF0ZSh0cmFuc2xvY2F0aW9uX2lkID0gYXMuZmFjdG9yKHRyYW5zbG9jYXRpb25faWQpKSAlPiUgCiAgaW5uZXJfam9pbihkMSwgYnkgPSBjKCJzaXRlX2lkIiwgInllYXIiKSkgJT4lIAogIHJlbG9jYXRlKHRyYW5zbG9jYXRpb25faWQsIC5hZnRlciA9IHNpdGVfaWQpCmBgYAoKIyMgTnVsbCB2YWx1ZXMgLSBjaGVjawoKYGBge3IgY2hlY2stbnVsbH0KZDEgJT4lIHN1bW1hcml6ZShhY3Jvc3MoZXZlcnl0aGluZygpLCB+IHN1bShpcy5uYSguKSkpKQpgYGAKCiMjIE1vZGVsIHVzaW5nIG5vbi1zdGFuZGFyZGl6ZWQgcHJlZGljdG9ycwojIyMgTW9kZWwgc3BlY2lmaWNhdGlvbgpgYGB7cn0KbTFhIDwtIHN0YW5fZ2xtKAogIHN1cnYgfiBsZW5ndGggKyBzbm93X3QgKyBzbm93X3QxICsgZGF5ICsgYmRsb2FkX2wgKyBlbGV2YXRpb24gKyBmaXJzdCArIG1hbGUgKyBkb25vciwKICBkYXRhID0gZDEsCiAgZmFtaWx5ID0gYmlub21pYWwsCiAgY2hhaW5zID0gNCwgCiAgaXRlciA9IDUwMDAqMiwKICBjb3JlcyA9IDQsCiAgc2VlZCA9IDg0NzM1KQpgYGAKCiMjIyBQcmlvcnMgJiBjaGFpbiBkaWFnbm9zdGljcwpgYGB7cn0KcHJpb3Jfc3VtbWFyeShtMWEpCgptY21jX3RyYWNlKG0xYSwgc2l6ZSA9IDAuMSkKbWNtY19kZW5zX292ZXJsYXkobTFhKQoKc3VtbWFyeShtMWEpCmBgYAotIG1jbWMgZGlhZ25vc3RpY3MgYWxsIGxvb2sgZ29vZC4KCiMjIyBMaXN0IG1vZGVsIHZhcmlhYmxlIG5hbWVzCmBgYHtyfQpnZXRfdmFyaWFibGVzKG0xYSkKYGBgCgojIyMgUG9zdGVyaW9yIHByZWRpY3RpdmUgY2hlY2sKYGBge3J9CnBwX2NoZWNrKG0xYSwgcGxvdGZ1biA9ICJzdGF0IikgKwogICAgeGxhYigiRnJvZyBzdXJ2aXZhbCByYXRlIikKYGBgCi0gUG9zdGVyaW9yIHByZWRpY3RpdmUgZGlzdHJpYnV0aW9uIGltcGxpZWQgYnkgcmVncmVzc2lvbiBtb2RlbCBjbG9zZWx5IG1hdGNoZXMgb3JpZ2luYWwgZGF0YSAoYmFzZWQgb24gbWVhbiBmcm9nIHN1cnZpdmFsKS4KCiMjIyBDYWxjdWxhdGUgbGVhdmUtb25lLW91dCBjcm9zcy12YWxpZGF0aW9uIChsb28pCmBgYHtyfQpsb29fbTFhIDwtIGxvbyhtMWEsIGNvcmVzID0gNCwgc2F2ZV9wc2lzID0gVFJVRSkKcHJpbnQobG9vX20xYSkKYGBgCgojIyMgUG9zdGVyaW9yIGFuYWx5c2lzCmBgYHtyfQp0aWR5KG0xYSwgY29uZi5pbnQgPSBUUlVFLCBjb25mLmxldmVsID0gMC45MCkKCm1jbWNfaW50ZXJ2YWxzKG0xYSwgcG9pbnRfZXN0ID0gIm1lZGlhbiIsIHByb2IgPSAwLjgwLCBwcm9iX291dGVyID0gMC45KSArCiAgZ2VvbV92bGluZSh4aW50ZXJjZXB0ID0gMCwgY29sb3IgPSAicmVkIiwgbGluZXR5cGUgPSAiZGFzaGVkIiwgbHdkID0gMC41KSArICAjIGFkZCAwLWxpbmUgdG8gYWNjZW50dWF0ZSBkZWZhdWx0IGxpbmUKICBnZ3RpdGxlKCJNb2RlbCB3aXRoIG5vbi1zdGFuZGFyZGl6ZWQgcHJlZGljdG9ycyBvbmx5IiwKICAgICAgICAgICJQb3N0ZXJpb3IgZGlzdHJpYnV0aW9ucyB3aXRoIG1lZGlhbnMgJiA5MCUgdW5jZXJ0YWludHkgaW50ZXJ2YWxzIikKCm1jbWNfYXJlYXMobTFhLCBhcmVhX21ldGhvZCA9ICJlcXVhbCBoZWlnaHQiLCBwcm9iID0gMC45MCkgKwogIGdlb21fdmxpbmUoeGludGVyY2VwdCA9IDAsIGNvbG9yID0gInJlZCIsIGxpbmV0eXBlID0gImRhc2hlZCIsIGx3ZCA9IDAuNSkgKyAgCiAgZ2d0aXRsZSgiTW9kZWwgd2l0aCBub24tc3RhbmRhcmRpemVkIHByZWRpY3RvcnMgb25seSIsCiAgICAgICAgICAiUG9zdGVyaW9yIGRpc3RyaWJ1dGlvbnMgd2l0aCBtZWRpYW5zICYgOTAlIHVuY2VydGFpbnR5IGludGVydmFscyIpCgptZWFuKGJheWVzX1IyKG0xYSkpCmBgYAotIFVzaW5nIG5vbi1zdGFuZGFyZGl6ZWQgcHJlZGljdG9ycyBtYWtlcyBwbG90IGhhcmQgdG8gcmVhZCBkdWUgdG8gdmFzdGx5IGRpZmZlcmVudCBzY2FsZXMgb2YgcHJlZGljdG9ycy4gIAotIEJhc2VkIG9uIDkwJSB1bmNlcnRhaW50eSBpbnRlcnZhbHMsIGltcG9ydGFudCBwcmVkaWN0b3JzIGFyZSBzbm93X3QxLCBlbGV2YXRpb24sIGZpcnN0LCBkb25vcjcwNTY3LCBhbmQgZG9ub3I3Mjk5Ni4KCiMjIE1vZGVsIHVzaW5nIHN0YW5kYXJkaXplZCBwcmVkaWN0b3JzCiMjIyBNb2RlbCBzcGVjaWZpY2F0aW9uCmBgYHtyfQptMWIgPC0gc3Rhbl9nbG0oCiAgc3VydiB+IGxlbmd0aF96ICsgc25vd3RfeiArIHNub3d0MV96ICsgZGF5X3ogKyBiZGxvYWRfeiArIGVsZXZfeiArIGZpcnN0X2MgKyBtYWxlX2MgKyBkb25vciwKICBkYXRhID0gZDEsCiAgZmFtaWx5ID0gYmlub21pYWwsCiAgY2hhaW5zID0gNCwgCiAgaXRlciA9IDUwMDAqMiwgCiAgY29yZXMgPSA0LCAKICBzZWVkID0gODQ3MzUpCmBgYAoKIyMjIFByaW9ycywgY2hhaW4gZGlhZ25vc3RpY3MsIHBvc3RlcmlvciBwcmVkaWN0aXZlIGNoZWNrCmBgYHtyfQpwcmlvcl9zdW1tYXJ5KG0xYikKCm1jbWNfdHJhY2UobTFiLCBzaXplID0gMC4xKQptY21jX2RlbnNfb3ZlcmxheShtMWEpCgpzdW1tYXJ5KG0xYikKYGBgCi0gbWNtYyBkaWFnbm9zdGljcyBhbGwgbG9vayBnb29kLgoKIyMjIExpc3QgbW9kZWwgdmFyaWFibGUgbmFtZXMKYGBge3J9CmdldF92YXJpYWJsZXMobTFiKQpgYGAKCiMjIyBQb3N0ZXJpb3IgcHJlZGljdGl2ZSBjaGVjawpgYGB7cn0KcHBfY2hlY2sobTFiLCBwbG90ZnVuID0gInN0YXQiKSArCiAgICB4bGFiKCJGcm9nIHN1cnZpdmFsIHJhdGUiKQpgYGAKLSBQb3N0ZXJpb3IgcHJlZGljdGl2ZSBkaXN0cmlidXRpb24gaW1wbGllZCBieSByZWdyZXNzaW9uIG1vZGVsIGNsb3NlbHkgbWF0Y2hlcyBvcmlnaW5hbCBkYXRhIChiYXNlZCBvbiBtZWFuIGZyb2cgc3Vydml2YWwpLgoKIyMjIENhbGN1bGF0ZSBsb28KYGBge3J9Cmxvb19tMWIgPC0gbG9vKG0xYiwgY29yZXMgPSA0LCBzYXZlX3BzaXMgPSBUUlVFKQpwcmludChsb29fbTFiKQpgYGAKCiMjIyBQb3N0ZXJpb3IgYW5hbHlzaXMKYGBge3J9CnRpZHkobTFiLCBjb25mLmludCA9IFRSVUUsIGNvbmYubGV2ZWwgPSAwLjkwKQoKbWNtY19hcmVhcyhtMWIsIGFyZWFfbWV0aG9kID0gImVxdWFsIGhlaWdodCIsIHByb2IgPSAwLjkwKSArCiAgZ2VvbV92bGluZSh4aW50ZXJjZXB0ID0gMCwgY29sb3IgPSAicmVkIiwgbGluZXR5cGUgPSAiZGFzaGVkIiwgbHdkID0gMC41KSArICAKICBnZ3RpdGxlKCJNb2RlbCB3aXRoIHN0YW5kYXJkaXplZCBwcmVkaWN0b3JzIG9ubHkiLAogICAgICAgICAgIlBvc3RlcmlvciBkaXN0cmlidXRpb25zIHdpdGggbWVkaWFucyAmIDkwJSB1bmNlcnRhaW50eSBpbnRlcnZhbHMiKQoKbWVhbihiYXllc19SMihtMWIpKQpgYGAKLSBBcyBleHBlY3RlZCwgdXNpbmcgc3RhbmRhcmRpemVkIHByZWRpY3RvcnMgYWxsb3dzIGVhc3kgY29tcGFyaXNvbiBhY3Jvc3MgcHJlZGljdG9ycy4gSW1wb3J0YW50IHByZWRpY3RvcnMgYXJlIHVuY2hhbmdlZC4gCi0gQmFzZWQgb24gOTAlIHVuY2VydGFpbnR5IGludGVydmFscywgaW1wb3J0YW50IHByZWRpY3RvcnMgYXJlIHNub3dfdDEsIGVsZXZhdGlvbiwgZmlyc3QsIGRvbm9yNzA1NjcsIGFuZCBkb25vcjcyOTk2LgoKIyMgQ29tcGFyZSBtb2RlbHMKYGBge3J9Cmxvb19jb21wYXJlKGxvb19tMWEsIGxvb19tMWIpCmBgYAotIEFzIGV4cGVjdGVkLCBtb2RlbHMgd2l0aCBub24tc3RhbmRhcmRpemVkIGFuZCBzdGFuZGFyZGl6ZWQgcHJlZGljdG9yIHZhcmlhYmxlcyBoYXZlIGlkZW50aWNhbCBwcmVkaWN0aXZlIGFjY3VyYWN5LiAgCgojIyBQb3NzaWJsZSBncm91cGluZyBzdHJ1Y3R1cmVzIHRvIGFkZCB0byBtb2RlbApUaGUgZGF0YXNldCBoYXMgdGhyZWUgZ3JvdXBzOiBzaXRlX2lkLCB0cmFuc2xvY2F0aW9uX2lkLCBhbmQgZG9ub3Igc2l0ZS4gRG9ub3Igc2l0ZSBpcyBub3Qgc3VpdGFibGUgYXMgYSBncm91cGluZyB2YXJpYWJsZSBiZWNhdXNlIHRoZXJlIGFyZSBvbmx5IDMgbGV2ZWxzIChpLmUuLCBzaXRlIGlkcyksIG5vdCBlbm91Z2ggdG8gbWFrZSB0aGlzIHVzZWZ1bC4gU2l0ZV9pZCBhbmQgdHJhbnNsb2NhdGlvbl9pZCBjb250YWluIDEyIGFuZCAyNCBsZXZlbHMsIHJlc3BlY3RpdmVseSwgbWFraW5nIHRoZW0gc3VpdGFibGUgYXMgZ3JvdXBpbmcgdmFyaWFibGVzLiBJbmNsdWRpbmcgYm90aCBzaXRlX2lkIGFuZCB0cmFuc2xvY2F0aW9uX2lkIGFzIG5lc3RlZCBncm91cGluZyB2YXJpYWJsZXMgKGkuZS4sIHRyYW5zbG9jYXRpb25zIG5lc3RlZCB3aXRoaW4gc2l0ZSBpZHMpIGFsc28gc2VlbXMganVzdGlmaWVkLiBIYXZpbmcgbW9yZSBsZXZlbHMgd2l0aGluIGEgZ3JvdXAgd291bGQgc2VlbSB0byBpbmNyZWFzZSB0aGUgcG93ZXIgdG8gZGV0ZWN0IGVmZmVjdHMgb2YgcHJlZGljdG9ycy4gV2l0aCBmZXdlciBsZXZlbHMvbGVzcyB2YXJpYWJpbGl0eSwgZWZmZWN0IG1heSBiZSBhY2NvdW50ZWQgZm9yIHByaW1hcmlseSBieSBncm91cGluZyB2YXJpYWJsZSBpdHNlbGYgaW5zdGVhZCBvZiBwcmVkaWN0b3IuIFdpdGggb25seSAxLTQgdHJhbnNsb2NhdGlvbnMgcGVyIHNpdGUgYW5kIG9ubHkgMTIgc2l0ZSBpZHMgdG90YWwsIGl0IGlzIHBvc3NpYmxlIHRoYXQgbW9kZWwgd2l0aCBzaXRlX2lkIGFzIHNvbGUgZ3JvdXBpbmcgdmFyaWFibGUgYW5kIG1vZGVsIHdpdGggc2l0ZV9pZCBhbmQgdHJhbnNsb2NhdGlvbl9pZCBhcyBuZXN0ZWQgZ3JvdXBpbmcgdmFyaWFibGVzIHdpbGwgYXR0cmlidXRlIG1vcmUgb2YgdGhlIHZhcmlhdGlvbiBpbiBmcm9nIHN1cnZpdmFsIHRvIHRoZSBncm91cGluZyB2YXJpYWJsZXMgYW5kIGxlc3MgdG8gdGhlIHByZWRpY3RvcnMuIEEgcHJpb3JpLCB0aGF0IHdvdWxkIHN1Z2dlc3QgdGhhdCB0aGUgbW9kZWwgd2l0aCB0cmFuc2xvY2F0aW9uX2lkIGFzIHRoZSBncm91cGluZyB2YXJpYWJsZSBtaWdodCBiZSB0aGUgbW9zdCB1c2VmdWwgbW9kZWwgc3RydWN0dXJlIChhc3N1bWluZyB0aGF0IGFsbCBtb2RlbHMgaGF2ZSBzaW1pbGFyIHBvc3RlcmlvciBwcmVkaWN0aXZlIGFjY3VyYWN5KS4gTGFja2luZyBhbnkgc3Ryb25nIHJhdGlvbmFsZSBmb3IgIG9uZSBtb2RlbCBzdHJ1Y3R1cmUgb3ZlciBhbm90aGVyLCBidWlsdCAzIG1vZGVscyBhbmQgY29tcGFyZWQgdGhlaXIgcG9zdGVyaW9yIHByZWRpY3RpdmUgYWNjdXJhY3kgdXNpbmcgbG9vLiAKCiMjIE1vZGVsIHVzaW5nIHN0YW5kYXJkaXplZCBwcmVkaWN0b3JzIGFuZCBzaXRlX2lkIGdyb3VwaW5nIHZhcmlhYmxlCiMjIyBNb2RlbCBzcGVjaWZpY2F0aW9uCmBgYHtyfQptMWMgPC0gc3Rhbl9nbG1lcigKICBzdXJ2IH4gbGVuZ3RoX3ogKyBzbm93dF96ICsgc25vd3QxX3ogKyBkYXlfeiArIGJkbG9hZF96ICsgZWxldl96ICsgZmlyc3RfYyArIG1hbGVfYyArIGRvbm9yICsgKDEgfCBzaXRlX2lkKSwKICBkYXRhID0gZDEsCiAgZmFtaWx5ID0gYmlub21pYWwsCiAgY2hhaW5zID0gNCwgCiAgaXRlciA9IDUwMDAqMiwgCiAgY29yZXMgPSA0LAogIHNlZWQgPSA4NDczNSkKYGBgCgojIyMgUHJpb3JzICYgY2hhaW4gZGlhZ25vc3RpY3MKYGBge3J9CnByaW9yX3N1bW1hcnkobTFjKQoKYjEgPC0gYygiKEludGVyY2VwdCkiLCAibGVuZ3RoX3oiLCAic25vd3RfeiIsICJzbm93dDFfeiIsICJkYXlfeiIsICJiZGxvYWRfeiIsICJlbGV2X3oiLCAiZmlyc3RfYzAuMzUzIiwgIm1hbGVfYzAuNTQ4IiwgImRvbm9yNzA1NjciLCAiZG9ub3I3Mjk5NiIpCgptY21jX3RyYWNlKG0xYywgc2l6ZSA9IDAuMSwgcGFycyA9IGIxKQptY21jX2RlbnNfb3ZlcmxheShtMWMsIHBhcnMgPSBiMSkKCnN1bW1hcnkobTFjKQpgYGAKCiMjIyBQb3N0ZXJpb3IgcHJlZGljdGl2ZSBjaGVjawpgYGB7cn0KcHBfY2hlY2sobTFjLCBwbG90ZnVuID0gInN0YXQiKSArCiAgICB4bGFiKCJGcm9nIHN1cnZpdmFsIHJhdGUiKQpgYGAKCiMjIyBDYWxjdWxhdGUgbG9vCmBgYHtyfQpsb29fbTFjIDwtIGxvbyhtMWMsIGNvcmVzID0gNCwgc2F2ZV9wc2lzID0gVFJVRSkKcHJpbnQobG9vX20xYykKYGBgCgojIyMgUG9zdGVyaW9yIGFuYWx5c2lzCmBgYHtyfQp0aWR5KG0xYywgY29uZi5pbnQgPSBUUlVFLCBjb25mLmxldmVsID0gMC45MCkKCm1jbWNfYXJlYXMobTFjLCBhcmVhX21ldGhvZCA9ICJlcXVhbCBoZWlnaHQiLCBwcm9iID0gMC45MCwgcGFycyA9IGIxKSArCiAgZ2VvbV92bGluZSh4aW50ZXJjZXB0ID0gMCwgY29sb3IgPSAicmVkIiwgbGluZXR5cGUgPSAiZGFzaGVkIiwgbHdkID0gMC41KSArICAKICBnZ3RpdGxlKCJNb2RlbCB3aXRoIHN0YW5kYXJkaXplZCBwcmVkaWN0b3JzICYgc2l0ZV9pZCBhcyBncm91cGluZyB2YXJpYWJsZSIsCiAgICAgICAgICAiUG9zdGVyaW9yIGRpc3RyaWJ1dGlvbnMgd2l0aCBtZWRpYW5zICYgOTAlIHVuY2VydGFpbnR5IGludGVydmFscyIpCgptZWFuKGJheWVzX1IyKG0xYykpCgojIEZyb2cgc3Vydml2YWwgYnkgc2l0ZV9pZAptY21jX2FyZWFzX3JpZGdlcyhtMWMsIHJlZ2V4X3BhcnMgPSAic2l0ZV9pZCIpICsKICBnZW9tX3ZsaW5lKHhpbnRlcmNlcHQgPSAwLCBjb2xvciA9ICJyZWQiLCBsaW5ldHlwZSA9ICJkYXNoZWQiLCBsd2QgPSAwLjUpICsgIAogIGdndGl0bGUoIkZyb2cgc3Vydml2YWwgZm9yIGVhY2ggcmVjaXBpZW50IHNpdGVfaWQiKQpgYGAKLSBCYXNlZCBvbiA5MCUgdW5jZXJ0YWludHkgaW50ZXJ2YWxzLCBpbXBvcnRhbnQgcHJlZGljdG9ycyByZWR1Y2VkIHRvIGVsZXZhdGlvbiBhbmQgbWFsZS4gCgojIyBNb2RlbCB1c2luZyBzdGFuZGFyZGl6ZWQgcHJlZGljdG9ycyBhbmQgdHJhbnNsb2NhdGlvbl9pZCBncm91cGluZyB2YXJpYWJsZSkKIyMjIE1vZGVsIHNwZWNpZmljYXRpb24KYGBge3J9Cm0xZCA8LSBzdGFuX2dsbWVyKAogIHN1cnYgfiBsZW5ndGhfeiArIHNub3d0X3ogKyBzbm93dDFfeiArIGRheV96ICsgYmRsb2FkX3ogKyBlbGV2X3ogKyBmaXJzdF9jICsgbWFsZV9jICsgZG9ub3IgKyAoMSB8IHRyYW5zbG9jYXRpb25faWQpLAogIGRhdGEgPSBkMSwKICBmYW1pbHkgPSBiaW5vbWlhbCwKICBjaGFpbnMgPSA0LCAKICBpdGVyID0gNTAwMCoyLCAKICBjb3JlcyA9IDQsCiAgc2VlZCA9IDg0NzM1KQpgYGAKCiMjIyBQcmlvcnMgJiBjaGFpbiBkaWFnbm9zdGljcwpgYGB7cn0KcHJpb3Jfc3VtbWFyeShtMWQpCgpiMSA8LSBjKCIoSW50ZXJjZXB0KSIsICJsZW5ndGhfeiIsICJzbm93dF96IiwgInNub3d0MV96IiwgImRheV96IiwgImJkbG9hZF96IiwgImVsZXZfeiIsICJmaXJzdF9jMC4zNTMiLCAibWFsZV9jMC41NDgiLCAiZG9ub3I3MDU2NyIsICJkb25vcjcyOTk2IikKCm1jbWNfdHJhY2UobTFkLCBzaXplID0gMC4xLCBwYXJzID0gYjEpCm1jbWNfZGVuc19vdmVybGF5KG0xZCwgcGFycyA9IGIxKQoKc3VtbWFyeShtMWQpCmBgYAoKIyMjIFBvc3RlcmlvciBwcmVkaWN0aXZlIGNoZWNrCmBgYHtyfQpwcF9jaGVjayhtMWQsIHBsb3RmdW4gPSAic3RhdCIpICsKICAgIHhsYWIoIkZyb2cgc3Vydml2YWwgcmF0ZSIpCmBgYAoKIyMjIENhbGN1bGF0ZSBsb28KYGBge3J9Cmxvb19tMWQgPC0gbG9vKG0xZCwgY29yZXMgPSA0LCBzYXZlX3BzaXMgPSBUUlVFKQpwcmludChsb29fbTFkKQpgYGAKCiMjIyBQb3N0ZXJpb3IgYW5hbHlzaXMKYGBge3J9CnRpZHkobTFkLCBjb25mLmludCA9IFRSVUUsIGNvbmYubGV2ZWwgPSAwLjkwKQoKbWNtY19hcmVhcyhtMWQsIGFyZWFfbWV0aG9kID0gImVxdWFsIGhlaWdodCIsIHByb2IgPSAwLjkwLCBwYXJzID0gYjEpICsKICBnZW9tX3ZsaW5lKHhpbnRlcmNlcHQgPSAwLCBjb2xvciA9ICJyZWQiLCBsaW5ldHlwZSA9ICJkYXNoZWQiLCBsd2QgPSAwLjUpICsgIAogIGdndGl0bGUoIk1vZGVsIHdpdGggc3RhbmRhcmRpemVkIHByZWRpY3RvcnMgJiB0cmFuc2xvY2F0aW9uX2lkIGFzIGdyb3VwaW5nIHZhcmlhYmxlIiwKICAgICAgICAgICJQb3N0ZXJpb3IgZGlzdHJpYnV0aW9ucyB3aXRoIG1lZGlhbnMgJiA5MCUgdW5jZXJ0YWludHkgaW50ZXJ2YWxzIikKCm1lYW4oYmF5ZXNfUjIobTFkKSkKCiMgRnJvZyBzdXJ2aXZhbCBieSB0cmFuc2xvY2F0aW9uX2lkCm1jbWNfYXJlYXNfcmlkZ2VzKG0xZCwgcmVnZXhfcGFycyA9ICJ0cmFuc2xvY2F0aW9uX2lkIikgKwogIGdlb21fdmxpbmUoeGludGVyY2VwdCA9IDAsIGNvbG9yID0gInJlZCIsIGxpbmV0eXBlID0gImRhc2hlZCIsIGx3ZCA9IDAuNSkgKyAgCiAgZ2d0aXRsZSgiRnJvZyBzdXJ2aXZhbCBmb3IgZWFjaCB0cmFuc2xvY2F0aW9uX2lkIikKYGBgCi0gQmFzZWQgb24gbG9vLCBtb2RlbCB3aXRoIHRyYW5zbG9jYXRpb25faWQgYXMgZ3JvdXBpbmcgdmFyaWFibGUgcHJvdmlkZXMgc2xpZ2h0bHkgYmV0dGVyIGZpdCB0aGFuIG1vZGVsIHdpdGggc2l0ZV9pZCBhcyBncm91cGluZyB2YXJpYWJsZS4gCi0gQmFzZWQgb24gOTAlIHVuY2VydGFpbnR5IGludGVydmFscywgaW1wb3J0YW50IHByZWRpY3RvcnMgYXJlIHNpbWlsYXIgdG8gdGhvc2UgaW4gbW9kZWwgd2l0aG91dCBncm91cGluZyB2YXJpYWJsZTogc25vd3QxX3osIGVsZXZfeiwgbWFsZV9jMC41NDgsIGRvbm9yNzI5OTYuIAoKIyMgTW9kZWwgdXNpbmcgc3RhbmRhcmRpemVkIHByZWRpY3RvcnMgYW5kIDIgZ3JvdXBpbmcgdmFyaWFibGVzOiB0cmFuc2xvY2F0aW9uX2lkIG5lc3RlZCB3aXRoaW4gc2l0ZV9pZCkKIyMjIE1vZGVsIHNwZWNpZmljYXRpb24KYGBge3J9Cm0xZSA8LSBzdGFuX2dsbWVyKAogIHN1cnYgfiBsZW5ndGhfeiArIHNub3d0X3ogKyBzbm93dDFfeiArIGRheV96ICsgYmRsb2FkX3ogKyBlbGV2X3ogKyBmaXJzdF9jICsgbWFsZV9jICsgZG9ub3IgKyAoMSB8IHRyYW5zbG9jYXRpb25faWQpICsgKDEgfCBzaXRlX2lkKSwKICBkYXRhID0gZDEsCiAgZmFtaWx5ID0gYmlub21pYWwsCiAgY2hhaW5zID0gNCwgCiAgaXRlciA9IDUwMDAqMiwgCiAgY29yZXMgPSA0LAogIHNlZWQgPSA4NDczNSkKYGBgCgojIyMgUHJpb3JzICYgY2hhaW4gZGlhZ25vc3RpY3MKYGBge3J9CnByaW9yX3N1bW1hcnkobTFlKQoKYjEgPC0gYygiKEludGVyY2VwdCkiLCAibGVuZ3RoX3oiLCAic25vd3RfeiIsICJzbm93dDFfeiIsICJkYXlfeiIsICJiZGxvYWRfeiIsICJlbGV2X3oiLCAiZmlyc3RfYzAuMzUzIiwgIm1hbGVfYzAuNTQ4IiwgImRvbm9yNzA1NjciLCAiZG9ub3I3Mjk5NiIpCgptY21jX3RyYWNlKG0xZSwgc2l6ZSA9IDAuMSwgcGFycyA9IGIxKQptY21jX2RlbnNfb3ZlcmxheShtMWUsIHBhcnMgPSBiMSkKCnN1bW1hcnkobTFlKQpgYGAKCiMjIyBQb3N0ZXJpb3IgcHJlZGljdGl2ZSBjaGVjawpgYGB7cn0KcHBfY2hlY2sobTFlLCBwbG90ZnVuID0gInN0YXQiKSArCiAgICB4bGFiKCJGcm9nIHN1cnZpdmFsIHJhdGUiKQpgYGAKCiMjIyBDYWxjdWxhdGUgbG9vCmBgYHtyfQpsb29fbTFlIDwtIGxvbyhtMWUsIGNvcmVzID0gNCwgc2F2ZV9wc2lzID0gVFJVRSkKcHJpbnQobG9vX20xZSkKYGBgCgojIyMgUG9zdGVyaW9yIGFuYWx5c2lzCmBgYHtyfQp0aWR5KG0xZSwgY29uZi5pbnQgPSBUUlVFLCBjb25mLmxldmVsID0gMC45MCkKCm1jbWNfYXJlYXMobTFlLCBhcmVhX21ldGhvZCA9ICJlcXVhbCBoZWlnaHQiLCBwcm9iID0gMC45MCwgcGFycyA9IGIxKSArCiAgZ2VvbV92bGluZSh4aW50ZXJjZXB0ID0gMCwgY29sb3IgPSAicmVkIiwgbGluZXR5cGUgPSAiZGFzaGVkIiwgbHdkID0gMC41KSArICAKICBnZ3RpdGxlKCJNb2RlbCB3aXRoIHN0YW5kYXJkaXplZCBwcmVkaWN0b3JzICYgdHJhbnNsb2NhdGlvbl9pZCBhcyBncm91cGluZyB2YXJpYWJsZSIsCiAgICAgICAgICAiUG9zdGVyaW9yIGRpc3RyaWJ1dGlvbnMgd2l0aCBtZWRpYW5zICYgOTAlIHVuY2VydGFpbnR5IGludGVydmFscyIpCgptZWFuKGJheWVzX1IyKG0xZSkpCgojIEZyb2cgc3Vydml2YWwgYnkgdHJhbnNsb2NhdGlvbl9pZAptY21jX2FyZWFzX3JpZGdlcyhtMWUsIHJlZ2V4X3BhcnMgPSAoInRyYW5zbG9jYXRpb25faWQiKSkgKwogIGdlb21fdmxpbmUoeGludGVyY2VwdCA9IDAsIGNvbG9yID0gInJlZCIsIGxpbmV0eXBlID0gImRhc2hlZCIsIGx3ZCA9IDAuNSkgKyAgCiAgZ2d0aXRsZSgiRnJvZyBzdXJ2aXZhbCBmb3IgZWFjaCB0cmFuc2xvY2F0aW9uX2lkIikKCiMgRnJvZyBzdXJ2aXZhbCBieSBzaXRlX2lkCm1jbWNfYXJlYXNfcmlkZ2VzKG0xZSwgcmVnZXhfcGFycyA9ICgic2l0ZV9pZCIpKSArCiAgZ2VvbV92bGluZSh4aW50ZXJjZXB0ID0gMCwgY29sb3IgPSAicmVkIiwgbGluZXR5cGUgPSAiZGFzaGVkIiwgbHdkID0gMC41KSArICAKICBnZ3RpdGxlKCJGcm9nIHN1cnZpdmFsIGZvciBlYWNoIHRyYW5zbG9jYXRpb25faWQiKQpgYGAKCiMjIENvbXBhcmUgNCBtb2RlbHMgd2l0aCBkaWZmZXJlbnQgZ3JvdXBpbmcgc3RydWN0dXJlcwpgYGB7cn0KbG9vX2NvbXBhcmUobG9vX20xYiwgbG9vX20xYywgbG9vX20xZCwgbG9vX20xZSkKYGBgCi0gQmFzZWQgb24gZXhwZWN0ZWQgbG9nLXByZWRpY3RpdmUgZGVuc2l0aWVzIChFTFBEKSwgMyBtb2RlbHMgd2l0aCBncm91cGluZyB2YXJpYWJsZXMgYWxsIGhhdmUgc3Vic3RhbnRpYWxseSBoaWdoZXIgcG9zdGVyaW9yIHByZWRpY3RpdmUgYWNjdXJhY3kgdGhhbiBtb2RlbCB3aXRob3V0IGFueSBncm91cGluZyB2YXJpYWJsZXMgKGVscGRfZGlmZiBvZiBtMWIgaXMgbW9yZSB0aGFuIDIgU0UgbG93ZXIgdGhhbiBlbHBkX2RpZmYgb2YgbTFjOiBlbHBkX2RpZmYgPSAtMTguMiAkXHBtJCAxMi42LCAyIFNFIGludGVydmFsIGRvZXMgbm90IG92ZXJsYXAgMDogLTUuNiAtLSAtMzAuOCkuIAotIDMgbW9kZWxzIHdpdGggMSBvciBtb3JlIGdyb3VwaW5nIHZhcmlhYmxlcyBhcmUgZXF1aXZhbGVudCBhbmQgaGF2ZSBzaW1pbGFyIHBvc3RlcmlvciBwcmVkaWN0aXZlIGFjY3VyYWN5IChlbHBkX2RpZmYgb2YgbTFlIGFuZCBtMWQgYXJlIHdpdGhpbiAyIFNFIG9mIGVscGRfZGlmZiBvZiBtMWMpLiAgCi0gbTFkIHNlZW1zIHRvIHByb3ZpZGUgbW9yZSBpbnNpZ2h0IGludG8gdGhlIGVmZmVjdCBvZiB0aGUgcHJlZGljdG9ycywgc28gY2hvc2UgdGhpcyBtb2RlbCBmb3IgYWRkaXRpb25hbCBhbmFseXNpcy4gCgojIyBBZGRpdGlvbmFsIGFuYWx5c2lzIG9mIG1vZGVsIG0xZAojIyMgUG9zdGVyaW9yIGNsYXNzaWZpY2F0aW9uIGFjY3VyYWN5CmBgYHtyfQojIFBvc3RlcmlvciBwcmVkaWN0aXZlIG1vZGVscyBmb3IgZWFjaCBmcm9nIGluIGRhdGFzZXQKc2V0LnNlZWQoODQ3MzUpCnN1cnZpdmFsX3ByZWRpY3RfMSA8LSBwb3N0ZXJpb3JfcHJlZGljdChtMWQsIG5ld2RhdGEgPSBkMSkKCiMgQ2xhc3NpZnkgZnJvZyBzdXJ2aXZhbCBiYXNlZCBvbiBtZWFuIHByZWRpY3RlZCBzdXJ2aXZhbCBwcm9iYWJpbGl0eQpzdXJ2aXZhbF9wcmVkaWN0XzIgPC0gZDEgJT4lIAogIG11dGF0ZShzdXJ2aXZhbF9wcm9iID0gY29sTWVhbnMoc3Vydml2YWxfcHJlZGljdF8xKSwKICAgICAgICAgc3Vydml2YWxfY2xhc3MgPSBhcy5udW1lcmljKHN1cnZpdmFsX3Byb2IgPj0gMC41KSkgJT4lIAogIHNlbGVjdChzaXRlX2lkLCB0cmFuc2xvY2F0aW9uX2lkLCB5ZWFyLCBzbm93dDFfeiwgZWxldl96LCBtYWxlX2MsIGRvbm9yLCBzdXJ2LCBzdXJ2aXZhbF9wcm9iLCBzdXJ2aXZhbF9jbGFzcykKCiMgR2VuZXJhdGUgY29uZnVzaW9uIG1hdHJpeCBhbmQgY2FsY3VsYXRlIGFjY3VyYWN5CnN1cnZpdmFsX3ByZWRpY3RfMiAlPiUgCiAgdGFieWwoc3Vydiwgc3Vydml2YWxfY2xhc3MpICU+JSAKICBhZG9ybl90b3RhbHMoYygicm93IiwgImNvbCIpKQoKKDQ3MSArIDEzOCkvNzY3ICMgb3ZlcmFsbF9hY2N1cmFjeSAKNDcxLzUzMyAjIHRydWUgbmVnYXRpdmUgcmF0ZQoxMzgvMjM0ICMgdHJ1ZSBwb3NpdGl2ZSByYXRlCgojIEdlbmVyYXRlIGNvbmZ1c2lvbiBtYXRyaXggdXNpbmcgZnVuY3Rpb24Kc2V0LnNlZWQoODQ3MzUpCmNsYXNzaWZpY2F0aW9uX3N1bW1hcnkobTFkLCBkYXRhID0gZDEsIGN1dG9mZiA9IDAuNSkKYGBgCi0gU2Vuc2l0aXZpdHkgb3IgInRydWUgcG9zaXRpdmUgcmF0ZSIgbWVhc3VyZXMgdGhlIHByb3BvcnRpb24gb2YgWSA9IDEgb2JzZXJ2YXRpb25zIHRoYXQgYXJlIGFjY3VyYXRlbHkgY2xhc3NpZmllZC4gIAotIFNwZWNpZmljaXR5IG9yICJ0cnVlIG5lZ2F0aXZlIHJhdGUiIG1lYXN1cmVzIHRoZSBwcm9wb3J0aW9uIG9mIFkgPSAwIG9ic2VydmF0aW9ucyB0aGF0IGFyZSBhY2N1cmF0ZWx5IGNsYXNzaWZpZWQgKGkuZS4sIHByb3BvcnRpb24gb2YgZnJvZ3MgY29ycmVjdGx5IGlkZW50aWZpZWQgYXMgbm90IHN1cnZpdmluZykuICAKCiMjIyBDcmVhdGUgcmVkdWNlZCBtb2RlbCBjb250YWluaW5nIG9ubHkgaW1wb3J0YW50IHByZWRpY3RvcnMgZm9yIHN1YnNlcXVlbnQgYW5hbHlzZXMKLSBSZWR1Y2VkIG1vZGVsIGlzIGhlbHBmdWwgYmVjYXVzZSBpdCByZWR1Y2VzIHRoZSBudW1iZXIgb2YgdmFsdWVzIHVzZWQgdG8gY2FsY3VsYXRlIHByZWRpY3RpdmUgbW9kZWxzLiAKYGBge3J9Cm0xZF9yZWR1Y2UgPC0gc3Rhbl9nbG1lcigKICBzdXJ2IH4gc25vd3QxX3ogKyBlbGV2X3ogKyBtYWxlX2MgKyBkb25vciArICgxIHwgdHJhbnNsb2NhdGlvbl9pZCksCiAgZGF0YSA9IGQxLAogIGZhbWlseSA9IGJpbm9taWFsLAogIGNoYWlucyA9IDQsIAogIGl0ZXIgPSA1MDAwKjIsIAogIGNvcmVzID0gNCwKICBzZWVkID0gODQ3MzUpCmBgYAoKIyMjIENvbmZpcm0gZml0IG9mIHJlZHVjZWQgbW9kZWwgaXMgY29tcGFyYWJsZSB0byBmdWxsIG1vZGVsCmBgYHtyfQpsb29fbTFkX3JlZHVjZSA8LSBsb28obTFkLCBjb3JlcyA9IDQsIHNhdmVfcHNpcyA9IFRSVUUpCnByaW50KGxvb19tMWQpCgpsb29fY29tcGFyZShsb29fbTFkLCBsb29fbTFkX3JlZHVjZSkKYGBgCgojIyMgU2ltdWxhdGUgYW5kIHBsb3QgcG9zdGVyaW9yIHByZWRpY3RpdmUgbW9kZWxzICYgb2JzZXJ2YXRpb25zCmBgYHtyfQptMWRfZ3JvdXBfbWVhbnMgPC0gZDEgJT4lIAogIGdyb3VwX2J5KHRyYW5zbG9jYXRpb25faWQpICU+JSAKICBzdW1tYXJpemUoY291bnQgPSBuKCksIHBzdXJ2ID0gbWVhbihzdXJ2KSwgZWxldl96ID0gbWVhbihlbGV2X3opLCBzbm93dDFfeiA9IG1lYW4oc25vd3QxX3opKSAlPiUgCiAgbXV0YXRlKAogICAgbWFsZV9jID0gYXMuZmFjdG9yKDAuNTQ4KSwgCiAgICBkb25vciA9IGNhc2Vfd2hlbigKICAgICAgc3RyX2RldGVjdCh0cmFuc2xvY2F0aW9uX2lkLCAiNzQ5NzYiKSB8IHN0cl9kZXRlY3QodHJhbnNsb2NhdGlvbl9pZCwgIjcwNTU2IikgfiAiNzA0NTkiLAogICAgICBzdHJfZGV0ZWN0KHRyYW5zbG9jYXRpb25faWQsICI3MDQxMyIpIH4gIjcwNTY3IiwKICAgICAgVFJVRSB+ICI3Mjk5NiIpKSAlPiUgCiAgYXJyYW5nZShwc3VydikgICMgdG8gYXJyYW5nZSBwbG90IHZhbHVlcy9sYWJlbHMgYnkgcHN1cnYKCnByZWRpY3Rpb25zX20xZCA8LSBwb3N0ZXJpb3JfcHJlZGljdChtMWRfcmVkdWNlLCBuZXdkYXRhID0gbTFkX2dyb3VwX21lYW5zKQoKcHBjX2ludGVydmFscygKICB5ID0gbTFkX2dyb3VwX21lYW5zJHBzdXJ2LAogIHlyZXAgPSBwcmVkaWN0aW9uc19tMWQsIAogIHByb2IgPSAwLjUsIAogIHByb2Jfb3V0ZXIgPSAwLjgpICsKZ2dwbG90Mjo6c2NhbGVfeF9jb250aW51b3VzKAogIGxhYmVscyA9IG0xZF9ncm91cF9tZWFucyR0cmFuc2xvY2F0aW9uX2lkLAogIGJyZWFrcyA9IDE6bnJvdyhtMWRfZ3JvdXBfbWVhbnMpKSArCiAgeGF4aXNfdGV4dChhbmdsZSA9IDkwLCB2anVzdCA9IDAuNSkgKwp4bGFiKCJUcmFuc2xvY2F0aW9uX2lkIikKYGBgCi0gVW5jZXJ0YWludHkgaW50ZXJ2YWxzIGFwcGVhciBpbmNvcnJlY3QgZm9yIGJpbm9taWFsIG1vZGVsLCBidXQgbm90IHN1cmUgaG93IHRvIGFkanVzdC4gIAotIEJhc2VkIG9uIG1lYW4gdmFsdWVzLCBtb2RlbCBkb2VzIGFuIGFkZXF1YXRlIGpvYiBvZiBwcmVkaWN0aW5nIHN1cnZpdmFsIGJhc2VkIG9uIGdyb3VwLWxldmVsIGFuZCBpbmRpdmlkdWFsLWxldmVsIGNoYXJhY3RlcmlzdGljcy4gIAotIFdoYXQgZXhwbGFpbnMgdGhlIGxvdyBwcmVkaWN0ZWQgc3Vydml2YWwgYXQgNzAzNzAgKGJvdGggdHJhbnNsb2NhdGlvbnMpLCBhbmQgaGlnaCB2ZXJzdXMgbG93IHByZWRpY3RlZCBzdXJ2aXZhbCBhdCA3MDU1MCBmb3IgMXN0IGFuZCAybmQgdHJhbnNsb2NhdGlvbnMsIHJlc3BlY3RpdmVseT8gQm90aCBzaXRlcyBzZWVtIHRvIGNvbnRhaW4gaGlnaCBxdWFsaXR5IGhhYml0YXQgKGkuZS4sIGhpZ2ggZWxldmF0aW9uKS4gIAotIFRoZSBwcmVkaWN0ZWQgc3Vydml2YWwgaXMgYmFzZWQgb24gYXNzdW1wdGlvbiB0aGF0IG9ubHkgbWFsZXMgYXJlIGluY2x1ZGVkLiBSdW4gYW5vdGhlciBzZXQgb2YgbW9kZWxzIGZvciBmZW1hbGVzPyAgCgojIyMgUGxvdCBwb3N0ZXJpb3IgcGxhdXNpYmxlIHJlbGF0aW9uc2hpcHMgYmV0d2VlbiBlYWNoIHByZWRpY3RvciBhbmQgdGhlIHByb2JhYmlsaXR5IG9mIGZyb2cgc3Vydml2YWwgKGNvbmRpdGlvbmFsIGVmZmVjdHMpCi0gQXMgZGVzY3JpYmVkIGluIGh0dHBzOi8vcmRyci5pby9jcmFuL2JybXMvbWFuL2NvbmRpdGlvbmFsX2VmZmVjdHMuYnJtc2ZpdC5odG1sLCAiV2hlbiBjcmVhdGluZyBjb25kaXRpb25hbF9lZmZlY3RzIGZvciBhIHBhcnRpY3VsYXIgcHJlZGljdG9yIChvciBpbnRlcmFjdGlvbiBvZiB0d28gcHJlZGljdG9ycyksIG9uZSBoYXMgdG8gY2hvb3NlIHRoZSB2YWx1ZXMgb2YgYWxsIG90aGVyIHByZWRpY3RvcnMgdG8gY29uZGl0aW9uIG9uLiBCeSBkZWZhdWx0LCB0aGUgbWVhbiBpcyB1c2VkIGZvciBjb250aW51b3VzIHZhcmlhYmxlcyBhbmQgdGhlIHJlZmVyZW5jZSBjYXRlZ29yeSBpcyB1c2VkIGZvciBmYWN0b3JzLCBidXQgeW91IG1heSBjaGFuZ2UgdGhlc2UgdmFsdWVzIHZpYSBhcmd1bWVudCBjb25kaXRpb25zLiBUaGlzIGFsc28gaGFzIGFuIGltcGxpY2F0aW9uIGZvciB0aGUgcG9pbnRzIGFyZ3VtZW50OiBJbiB0aGUgY3JlYXRlZCBwbG90cywgb25seSB0aG9zZSBwb2ludHMgd2lsbCBiZSBzaG93biB0aGF0IGNvcnJlc3BvbmQgdG8gdGhlIGZhY3RvciBsZXZlbHMgYWN0dWFsbHkgdXNlZCBpbiB0aGUgY29uZGl0aW9uaW5nLCBpbiBvcmRlciBub3QgdG8gY3JlYXRlIHRoZSBmYWxzZSBpbXByZXNzaW9uIG9mIGJhZCBtb2RlbCBmaXQsIHdoZXJlIGl0IGlzIGp1c3QgZHVlIHRvIGNvbmRpdGlvbmluZyBvbiBjZXJ0YWluIGZhY3RvciBsZXZlbHMuIiAgCi0gVXNpbmcgcnN0YW5hcm0sIEkgd2FzIHVuYWJsZSB0byBkbyB0aGlzIGFzIEkgZXhwZWN0ZWQgdG86IFN0YXJ0IHdpdGggZGF0YSA9IGQxIGFuZCB1c2luZyAibmV3ZGF0YSA9ICIgdG8gc3BlY2lmeSB0aGF0IHByZWRpY3Rpb25zIHNob3VsZCBiZSBiYXNlZCBvbiBkYXRhIGluIGQxX2VsZXZ6IChmb3IgZXhhbXBsZSksIGluIHdoaWNoIG5vbi1wbG90dGVkIGNvbnRpbnVvdXMgcHJlZGljdG9ycyBhcmUgc2V0IHRvIHRoZSBtZWFuIHZhbHVlLiBOb3Qgc3VyZSBpZiBteSBhbHRlcm5hdGl2ZSBhcHByb2FjaCBpcyBjb3JyZWN0LiAgCi0gQXBwcm9hY2ggdXNlZDogSW5zdGVhZCBvZiBzdGFydGluZyB0aGUgcm91dGluZSB3aXRoIGQxLCBJIHN0YXJ0ZWQgd2l0aCB0aGUgbmV3IGRhdGFzZXQgaW4gd2hpY2ggYWxsIHZhcmlhYmxlcyBleGNlcHQgdGhlIG9uZSBvZiBpbnRlcmVzdCB3ZXJlIHNldCB0byBtZWFuIHZhbHVlcyAoY29udGludW91cyB2YXJpYWJsZXMpIGFuZCB0aGUgcmVmZXJlbmNlIGxldmVsIChjYXRlZ29yaWNhbCB2YXJpYWJsZXMpLiBUaGlzIHByb2R1Y2VkIHBsb3RzL3Jlc3VsdHMgdGhhdCB3ZXJlIHZlcnkgc2ltaWxhciB0byB0aG9zZSBwcm9kdWNlZCBpbiBicm1zIHVzaW5nIHRoZSBjb25kaXRpb25hbF9lZmZlY3RzIGZ1bmN0aW9uLiAKCiMjIyMgRWZmZWN0IG9mIGVsZXZhdGlvbgpgYGB7cn0KZDFfZWxldnogPC0gZDEgJT4lIAogIG11dGF0ZShzbm93dDFfeiA9IG1lYW4oc25vd3QxX3opLCAKICAgICAgICAgbWFsZV9jID0gcmVjb2RlKG1hbGVfYywgIi0wLjQ1MiIgPSAiMC41NDgiKSwKICAgICAgICAgZG9ub3IgPSByZWNvZGUoZG9ub3IsICI3MDU2NyIgPSAiNzA0NTkiLCAiNzI5OTYiID0gIjcwNDU5IikpCgpkMV9lbGV2eiAlPiUgZGlzdGluY3Qoc25vd3QxX3osIG1hbGVfYywgZG9ub3IpICMgZW5zdXJlIHZhcmlhYmxlcyB3ZXJlIHNldCB0byBjb3JyZWN0IHZhbHVlcwogIApkMV9lbGV2eiAlPiUgCiAgc2VsZWN0KHNpdGVfaWQsIHRyYW5zbG9jYXRpb25faWQsIHNub3d0MV96LCBlbGV2X3osIG1hbGVfYywgZG9ub3IsIHN1cnYpICU+JSAKICBhZGRfZXByZWRfZHJhd3MobTFkX3JlZHVjZSwgbmRyYXdzID0gMTAwLCByZV9mb3JtdWxhID0gTkEpICU+JSAKICBnZ3Bsb3QoYWVzKHggPSBlbGV2X3osIHkgPSBzdXJ2KSkgKwogICAgZ2VvbV9saW5lKGFlcyh5ID0gLmVwcmVkLCBncm91cCA9IC5kcmF3KSwgY29sb3IgPSAiYmx1ZSIsIHNpemUgPSAwLjIsIGFscGhhID0gMC41KSArCiMgICAgZ2VvbV9wb2ludChhZXMoeSA9IC5lcHJlZCkpICsKICBsYWJzKHggPSAiRWxldmF0aW9uIG9mIHJlY2lwaWVudCBzaXRlIChzdGFuZGFyZGl6ZWQpIiwgeSA9ICJQcm9iYWJpbGl0eSBvZiAxLXllYXIgZnJvZyBzdXJ2aXZhbCIpCgoKIyBkMSAlPiUgICMgVGVzdCBvZiBhZGRfcHJlZGljdGVkX2RyYXdzIHVzZWQgZm9yIHNpbWlsYXIgcHVycG9zZSAtIHVuc3VjY2Vzc2Z1bAojICAgc2VsZWN0KHNpdGVfaWQsIHRyYW5zbG9jYXRpb25faWQsIHNub3d0MV96LCBlbGV2X3osIG1hbGVfYywgZG9ub3IsIHN1cnYpICU+JSAKIyAgIGFkZF9wcmVkaWN0ZWRfZHJhd3MobTFkX3JlZHVjZSwgbmRyYXdzID0gMTAwLCByZV9mb3JtdWxhID0gTkEsIG5ld2RhdGEgPSBkMV9lbGV2eikgJT4lIFZpZXcoKQojICAgZ2dwbG90KGFlcyh4ID0gZWxldl96LCB5ID0gc3VydikpICsKIyAgICAgZ2VvbV9saW5lKGFlcyh5ID0gLmVwcmVkLCBncm91cCA9IC5kcmF3KSwgY29sb3IgPSAiYmx1ZSIsIHNpemUgPSAwLjIsIGFscGhhID0gMC41KSArCiMgICAgIGdlb21fcG9pbnQoYWVzKHkgPSAuZXByZWQpKSArCiMgICBsYWJzKHggPSAiRWxldmF0aW9uIG9mIHJlY2lwaWVudCBzaXRlIChzdGFuZGFyZGl6ZWQpIiwgeSA9ICJQcm9iYWJpbGl0eSBvZiAxLXllYXIgZnJvZyBzdXJ2aXZhbCIpCmBgYAoKIyMjIyBFZmZlY3Qgb2Ygd2ludGVyIHNldmVyaXR5IGluIHllYXIgdDEKYGBge3J9CmQxX3Nub3d0MXogPC0gZDEgJT4lIAogIG11dGF0ZShlbGV2X3ogPSBtZWFuKGVsZXZfeiksIAogICAgICAgICBtYWxlX2MgPSByZWNvZGUobWFsZV9jLCAiLTAuNDUyIiA9ICIwLjU0OCIpLAogICAgICAgICBkb25vciA9IHJlY29kZShkb25vciwgIjcwNTY3IiA9ICI3MDQ1OSIsICI3Mjk5NiIgPSAiNzA0NTkiKSkKCmQxX3Nub3d0MXogJT4lIGRpc3RpbmN0KGVsZXZfeiwgbWFsZV9jLCBkb25vcikgIyBlbnN1cmUgdmFyaWFibGVzIHdlcmUgc2V0IHRvIGNvcnJlY3QgdmFsdWVzCiAgCmQxX3Nub3d0MXogJT4lIAogIHNlbGVjdChzaXRlX2lkLCB0cmFuc2xvY2F0aW9uX2lkLCBzbm93dDFfeiwgZWxldl96LCBtYWxlX2MsIGRvbm9yLCBzdXJ2KSAlPiUgCiAgYWRkX2VwcmVkX2RyYXdzKG0xZF9yZWR1Y2UsIG5kcmF3cyA9IDEwMCwgcmVfZm9ybXVsYSA9IE5BKSAlPiUgCiAgZ2dwbG90KGFlcyh4ID0gc25vd3QxX3osIHkgPSBzdXJ2KSkgKwogICAgZ2VvbV9saW5lKGFlcyh5ID0gLmVwcmVkLCBncm91cCA9IC5kcmF3KSwgY29sb3IgPSAiYmx1ZSIsIHNpemUgPSAwLjIsIGFscGhhID0gMC41KSArCiMgICAgZ2VvbV9wb2ludChhZXMoeSA9IC5lcHJlZCkpICsKICBsYWJzKHggPSAiV2ludGVyIHNldmVyaXR5IGluIHllYXIgZm9sbG93aW5nIHRyYW5zbG9jYXRpb24gKHN0YW5kYXJkaXplZCkiLCB5ID0gIlByb2JhYmlsaXR5IG9mIDEteWVhciBmcm9nIHN1cnZpdmFsIikgCmBgYAoKIyMjIyBFZmZlY3Qgb2Ygc2V4CmBgYHtyfQpkMV9tYWxlYyA8LSBkMSAlPiUgCiAgbXV0YXRlKGVsZXZfeiA9IG1lYW4oZWxldl96KSwgCiAgICAgICAgIHNub3d0MV96ID0gbWVhbihzbm93dDFfeiksCiAgICAgICAgIGRvbm9yID0gcmVjb2RlKGRvbm9yLCAiNzA1NjciID0gIjcwNDU5IiwgIjcyOTk2IiA9ICI3MDQ1OSIpKQoKZDFfbWFsZWMgJT4lIGRpc3RpbmN0KGVsZXZfeiwgc25vd3QxX3osIGRvbm9yKSAjIGVuc3VyZSB2YXJpYWJsZXMgd2VyZSBzZXQgdG8gY29ycmVjdCB2YWx1ZXMKICAKZDFfbWFsZWMgJT4lIAogIHNlbGVjdChzaXRlX2lkLCB0cmFuc2xvY2F0aW9uX2lkLCBzbm93dDFfeiwgZWxldl96LCBtYWxlX2MsIGRvbm9yLCBzdXJ2KSAlPiUgCiAgYWRkX2VwcmVkX2RyYXdzKG0xZF9yZWR1Y2UsIG5kcmF3cyA9IDEwMCwgcmVfZm9ybXVsYSA9IE5BKSAlPiUgCiAgZ2dwbG90KGFlcyh4ID0gbWFsZV9jLCB5ID0gc3VydikpICsKICAgIGdlb21fYm94cGxvdChhZXMoeSA9IC5lcHJlZCwgeCA9IG1hbGVfYykpICsKIyAgICBnZW9tX3BvaW50KGFlcyh5ID0gLmVwcmVkKSkgKwogICAgbGFicyh4ID0gIlNleCAoY2VudGVyZWQpIiwgeSA9ICJQcm9iYWJpbGl0eSBvZiAxLXllYXIgZnJvZyBzdXJ2aXZhbCIpCmBgYAoKIyMjIyBFZmZlY3Qgb2YgZG9ub3IKYGBge3J9CmQxX2Rvbm9yIDwtIGQxICU+JSAKICBtdXRhdGUoZWxldl96ID0gbWVhbihlbGV2X3opLCAKICAgICAgICAgc25vd3QxX3ogPSBtZWFuKHNub3d0MV96KSwKICAgICAgICAgbWFsZV9jID0gcmVjb2RlKG1hbGVfYywgIi0wLjQ1MiIgPSAiMC41NDgiKSkKCmQxX2Rvbm9yICU+JSBkaXN0aW5jdChlbGV2X3osIHNub3d0MV96LCBtYWxlX2MpICMgZW5zdXJlIHZhcmlhYmxlcyB3ZXJlIHNldCB0byBjb3JyZWN0IHZhbHVlcwogIApkMV9kb25vciAlPiUgCiAgc2VsZWN0KHNpdGVfaWQsIHRyYW5zbG9jYXRpb25faWQsIHNub3d0MV96LCBlbGV2X3osIG1hbGVfYywgZG9ub3IsIHN1cnYpICU+JSAKICBhZGRfZXByZWRfZHJhd3MobTFkX3JlZHVjZSwgbmRyYXdzID0gMTAwLCByZV9mb3JtdWxhID0gTkEpICU+JSAKICBnZ3Bsb3QoYWVzKHggPSBkb25vciwgeSA9IHN1cnYpKSArCiAgICBnZW9tX2JveHBsb3QoYWVzKHkgPSAuZXByZWQsIHggPSBkb25vcikpICsKIyAgICBnZW9tX3BvaW50KGFlcyh5ID0gLmVwcmVkKSkgKwogICAgbGFicyh4ID0gIkRvbm9yIHBvcHVsYXRpb24iLCB5ID0gIlByb2JhYmlsaXR5IG9mIDEteWVhciBmcm9nIHN1cnZpdmFsIikKYGBgCgojIyMjIEpvaW50IGVmZmVjdHMgb2YgZWxldmF0aW9uIGFuZCBzZXgKYGBge3J9CmQxX2VsZXZ6X21hbGVjIDwtIGQxICU+JSAKICBtdXRhdGUoc25vd3QxX3ogPSBtZWFuKHNub3d0MV96KSwgCiAgICAgICAgIGRvbm9yID0gcmVjb2RlKGRvbm9yLCAiNzA1NjciID0gIjcwNDU5IiwgIjcyOTk2IiA9ICI3MDQ1OSIpKQoKZDFfZWxldnpfbWFsZWMgJT4lIGRpc3RpbmN0KHNub3d0MV96LCBkb25vcikgIyBlbnN1cmUgdmFyaWFibGVzIHdlcmUgc2V0IHRvIGNvcnJlY3QgdmFsdWVzCiAgCnBhbCA8LSBjKCIjMGQwODg3IiwgIiNmODk1NDAiKQpkMV9lbGV2el9tYWxlYyAlPiUgCiAgc2VsZWN0KHNpdGVfaWQsIHRyYW5zbG9jYXRpb25faWQsIHNub3d0MV96LCBlbGV2X3osIG1hbGVfYywgZG9ub3IsIHN1cnYpICU+JSAKICBhZGRfZXByZWRfZHJhd3MobTFkX3JlZHVjZSwgbmRyYXdzID0gMTAwLCByZV9mb3JtdWxhID0gTkEpICU+JSAKICBnZ3Bsb3QoYWVzKHggPSBlbGV2X3osIHkgPSBzdXJ2LCBjb2xvciA9IG1hbGVfYykpICsKICAgIGdlb21fbGluZShhZXMoeSA9IC5lcHJlZCwgZ3JvdXAgPSBwYXN0ZShtYWxlX2MsIC5kcmF3KSksIHNpemUgPSAwLjQsIGFscGhhID0gMC41KSArCiMgICAgZ2VvbV9wb2ludChhZXMoeSA9IC5lcHJlZCkpICsKICAgIGxhYnMoeCA9ICJFbGV2YXRpb24gb2YgcmVjaXBpZW50IHNpdGUgKHN0YW5kYXJkaXplZCkiLCB5ID0gIlByb2JhYmlsaXR5IG9mIDEteWVhciBmcm9nIHN1cnZpdmFsIikgKwogICAgc2NhbGVfY29sb3JfbWFudWFsKHZhbHVlcyA9IHBhbCkgICsgCiAgICB0aGVtZV9jbGFzc2ljKCkKYGBgCgoKIyMgRm9yIGNvbXBhcmlzb24sIGZpdCBtb2RlbCB1c2luZyBicm1zIGFuZCBjb25kaXRpb25hbF9lZmZlY3RzIGZ1bmN0aW9uIHRvIHBsb3QgdmFyaWFibGUtc3BlY2lmaWMgZWZmZWN0cwojIyMgRml0IG1vZGVsCmBgYHtyfQptMWRfcmVkdWNlX2JybXMgPC0gYnJtKAogIHN1cnYgfiBzbm93dDFfeiArIGVsZXZfeiArIG1hbGVfYyArIGRvbm9yICsgKDEgfCB0cmFuc2xvY2F0aW9uX2lkKSwKICBkYXRhID0gZDEsCiAgZmFtaWx5ID0gYmVybm91bGxpKCksCiAgY2hhaW5zID0gNCwgCiAgaXRlciA9IDUwMDAqMiwgCiAgY29yZXMgPSA0LAogIHNlZWQgPSA4NDczNSkKYGBgCgojIyMgUGxvdCBjb25kaXRpb25hbCBlZmZlY3RzCmBgYHtyfQpwbG90KGNvbmRpdGlvbmFsX2VmZmVjdHMobTFkX3JlZHVjZV9icm1zLCBzcGFnaGV0dGkgPSBUUlVFLCBtZWFuID0gVFJVRSwgbmRyYXdzID0gMTAwKSwgYXNrID0gRkFMU0UpCmBgYAoKCiMjIFRPIERPCi0gVmFsaWRhdGUgcHJpb3JzLCB0dW5lIGFzIG5lY2Vzc2FyeS4gQWRkIGZpbmFsIHByaW9ycyB0byBtb2RlbCBzcGVjaWZpY2F0aW9uIChpbiBjYXNlIHRoZSBkZWZhdWx0cyBjaGFuZ2UpLiAgCi0gQ2hlY2sgYXBwcm9hY2ggdG8gcGxvdHRpbmcgY29uZGl0aW9uYWwgZWZmZWN0cy4gIAotIENyZWF0ZSBhbmFseXNpcyBpbiB3aGljaCBzdXJ2aXZhbCBmcm9tIGZpcnN0IHRyYW5zbG9jYXRpb24gdG8gYSBzaXRlIGlzIGEgcHJlZGljdG9yIG9mIHN1cnZpdmFsIGZyb20gc3Vic2VxdWVudCB0cmFuc2xvY2F0aW9ucy4gRm9yIHRoaXMgYW5hbHlzaXMsIHdpbGwgbmVlZCB0byBkcm9wIGFsbCBzaXRlcyBhdCB3aGljaCBvbmx5IGEgc2luZ2xlIHRyYW5zbG9jYXRpb24gd2FzIGNvbmR1Y3RlZC4gIAotIFRyeSB1bi1zdGFuZGFyZGl6aW5nIHByZWRpY3RvciB2YWx1ZXMgaW4gY29uZGl0aW9uYWwtZWZmZWN0cyBwbG90cyAoYW5kIG90aGVyIHBsb3RzPykgdG8gcHJvZHVjZSBtb3JlIGludHVpdGl2ZSB2aXN1YWxpemF0aW9ucy4gIAoKCgoKCgoKCg==